SCM

SCM Repository

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

Diff of /pkg/src/dgBCMatrix.c

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

revision 535, Wed Feb 9 08:51:08 2005 UTC revision 536, Thu Feb 10 08:15:12 2005 UTC
# Line 27  Line 27 
27      return ScalarLogical(1);      return ScalarLogical(1);
28  }  }
29    
 static int*  
 expand_column_pointers(int ncol, const int mp[], int mj[])  
 {  
     int j;  
     for (j = 0; j < ncol; j++) {  
         int j2 = mp[j+1], jj;  
         for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;  
     }  
     return mj;  
 }  
   
 static R_INLINE  
 int Tind(const int rowind[], const int colptr[], int i, int j)  
 {  
     int k, k2 = colptr[j + 1];  
     for (k = colptr[j]; k < k2; k++)  
         if (rowind[k] == i) return k;  
     error("row %d and column %d not defined in rowind and colptr",  
           i, j);  
     return -1;                  /* -Wall */  
 }  
   
30  /**  /**
31   * Perform one of the matrix operations   * Perform one of the matrix operations
32   *  C := alpha*op(A)*B + beta*C   *  C := alpha*op(A)*B + beta*C
# Line 198  Line 176 
176          for (j = 0; j < anc; j++) {          for (j = 0; j < anc; j++) {
177              int k, kk, k2 = Ap[j+1];              int k, kk, k2 = Ap[j+1];
178              for (k = Ap[j]; k < k2; k++) {              for (k = Ap[j]; k < k2; k++) {
179                  int ii = Ai[k], K = Tind(Ci, Cp, ii, ii);                  int ii = Ai[k], K = check_csc_index(Cp, Ci, ii, ii, 0);
180    
181                  if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, ii);                  if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, ii);
182                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
# Line 207  Line 185 
185                                       &one, Cx + K * csz, cdims);                                       &one, Cx + K * csz, cdims);
186                  for (kk = k+1; kk < k2; kk++) {                  for (kk = k+1; kk < k2; kk++) {
187                      int jj = Ai[kk];                      int jj = Ai[kk];
188                      K = (uplo == UPP) ? Tind(Ci, Cp, ii, jj) : Tind(Ci, Cp, jj, ii);                      K = (uplo == UPP) ? check_csc_index(Cp, Ci, ii, jj, 0) :
189                            check_csc_index(Cp, Ci, jj, ii, 0);
190    
191                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];
192                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
# Line 222  Line 201 
201    
202  /**  /**
203   * Create the LD^{T/2}D^{1/2}L' decomposition of the positive definite   * Create the LD^{T/2}D^{1/2}L' decomposition of the positive definite
204   * symmetric dgBCMatrix matrix A (upper triangle stored) in L and   * symmetric dgBCMatrix matrix A (upper triangle stored) in L and D^{1/2}.
205   * D^{1/2}.  The notation D^{1/2} denotes the upper Cholesky factor of   * D^{1/2} denotes the upper Cholesky factor of the positive definite positive
206   * the positive definite positive definite block diagonal matrix D.   * definite block diagonal matrix D.  The diagonal blocks are of size nci.
  * The diagonal blocks are of size nci.  
207   *   *
208   * @param A pointer to a dgBCMatrix object containing the upper   * @param A pointer to a dgBCMatrix object containing the upper
209   * triangle of a positive definite symmetric matrix.   * triangle of a positive definite symmetric matrix.
# Line 248  Line 226 
226          *Ap = INTEGER(ApP),          *Ap = INTEGER(ApP),
227          *Li = INTEGER(GET_SLOT(L, Matrix_iSym)),          *Li = INTEGER(GET_SLOT(L, Matrix_iSym)),
228          *Lp = INTEGER(GET_SLOT(L, Matrix_pSym)), nci = adims[0];          *Lp = INTEGER(GET_SLOT(L, Matrix_pSym)), nci = adims[0];
229        int ncisqr = nci * nci;
230      double *Lx = REAL(GET_SLOT(L, Matrix_xSym)),      double *Lx = REAL(GET_SLOT(L, Matrix_xSym)),
231          *Ax = REAL(AxP), *Dx = REAL(D), minus1 = -1., one = 1.;          *Ax = REAL(AxP), *Dx = REAL(D), minus1 = -1., one = 1.;
232    
# Line 257  Line 236 
236          if (Parent[j] >= 0) {diag = 0; break;}          if (Parent[j] >= 0) {diag = 0; break;}
237      }      }
238      if (diag) {      if (diag) {
         int ncisqr = nci * nci;  
239          Memcpy(Dx, Ax, ncisqr * n);          Memcpy(Dx, Ax, ncisqr * n);
240          for (j = 0; j < n; j++) { /* form D_i^{1/2} */          for (j = 0; j < n; j++) { /* form D_i^{1/2} */
241              F77_CALL(dpotrf)("U", &nci, Dx + j * ncisqr, &nci, &k);              F77_CALL(dpotrf)("U", &nci, Dx + j * ncisqr, &nci, &k);
242                /* Should we just return k instead of signaling an error? */
243                /* The caller should check for return value < n */
244              if (k) error("D[ , , %d], level %d, is not positive definite",              if (k) error("D[ , , %d], level %d, is not positive definite",
245                           j + 1);                           j + 1);
246          }          }
# Line 274  Line 254 
254          return n;          return n;
255      } else {               /* Copy of ldl_numeric from the LDL package      } else {               /* Copy of ldl_numeric from the LDL package
256                              * modified for blocked sparse matrices */                              * modified for blocked sparse matrices */
257          int i, k, p, p2, len, nci = adims[0], ncisqr = adims[0]*adims[0], top;          int i, k, p, p2, len, top;
258          int *Lnz = Calloc(n, int),          int *Lnz = Calloc(n, int),
259              *Pattern = Calloc(n, int),              *Pattern = Calloc(n, int),
260              *Flag = Calloc(n, int);              *Flag = Calloc(n, int);
# Line 317  Line 297 
297                                      &one, Y + Li[p]*ncisqr, &nci);                                      &one, Y + Li[p]*ncisqr, &nci);
298                  }                  }
299                  /* FIXME: Check that this is the correct order and transposition */                  /* FIXME: Check that this is the correct order and transposition */
300                  F77_CALL(dtrsm)("R", "U", "N", "N", &nci, &nci,                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,
301                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
302                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,
303                                  &one, Dx + k*ncisqr, &nci);                                  &one, Dx + k*ncisqr, &nci);
304                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,                  F77_CALL(dtrsm)("R", "U", "N", "N", &nci, &nci,
305                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
306                  Li[p] = k;      /* store L[k,i] in column form of L */                  Li[p] = k;      /* store L[k,i] in column form of L */
307                  Memcpy(Lx + p * ncisqr, Yi, ncisqr);                  Memcpy(Lx + p * ncisqr, Yi, ncisqr);
# Line 336  Line 316 
316          Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);          Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);
317          return n;       /* success, diagonal of D is all nonzero */          return n;       /* success, diagonal of D is all nonzero */
318      }      }
319      return -1;                  /* keep -Wall happy */      return -1;                  /* -Wall */
320  }  }
321    
322  /**  /**
# Line 511  Line 491 
491  }  }
492    
493  /**  /**
  * Expand a column of a compressed, sparse, column-oriented matrix.  
  *  
  * @param dest array to hold the result  
  * @param m number of rows in the matrix  
  * @param j index (0-based) of column to expand  
  * @param Ap array of column pointers  
  * @param Ai array of row indices  
  * @param Ax array of non-zero values  
  *  
  * @return dest  
  */  
 static  
 double *expand_column(double *dest, int m, int j,  
                       const int Ap[], const int Ai[], const double Ax[])  
 {  
     int k, k2 = Ap[j + 1];  
   
     for (k = 0; k < m; k++) dest[k] = 0.;  
     for (k = Ap[j]; k < k2; k++) dest[Ai[k]] = Ax[k];  
     return dest;  
 }  
   
 /**  
494   * Solve one of the systems op(A)*X = alpha*B or   * Solve one of the systems op(A)*X = alpha*B or
495   * X*op(A) = alpha*B where A dgBCMatrix triangular and B is dgBCMatrix.   * X*op(A) = alpha*B where A dgBCMatrix triangular and B is dgBCMatrix.
496   *   *
# Line 591  Line 548 
548              rhs = Calloc(ncbB, double);              rhs = Calloc(ncbB, double);
549              AZERO(Bx, nnz);     /* zero the result */              AZERO(Bx, nnz);     /* zero the result */
550              for (i = 0; i < nrbB; i++) {              for (i = 0; i < nrbB; i++) {
551                  R_ldl_lsolve(ncbB, expand_column(rhs, ncbB, i, BTp, BTi, BTx),                  R_ldl_lsolve(ncbB, expand_csc_column(rhs, ncbB, i, BTp, BTi, BTx),
552                               Ap, Ai, Ax);                               Ap, Ai, Ax);
553                  for (j = 0; j < ncbB; j++) { /* write non-zeros in sol'n into B */                  for (j = 0; j < ncbB; j++) { /* write non-zeros in sol'n into B */
554                      if (BTx[j]) Bx[Tind(Bi, Bp, j, i)] = BTx[j];                      if (BTx[j]) Bx[check_csc_index(Bp, Bi, j, i, 0)] = BTx[j];
555                  }                  }
556                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);
557              }              }
# Line 668  Line 625 
625              F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,              F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
626                                      &alpha, Ax + ia * ablk, adims,                                      &alpha, Ax + ia * ablk, adims,
627                                      Bx + ib * bblk, bdims, &one,                                      Bx + ib * bblk, bdims, &one,
628                                      Cx + Tind(Ci, Cp, Ai[ia], Bi[ib])*cblk, cdims);                              Cx + check_csc_index(Cp, Ci, Ai[ia], Bi[ib], 0) * cblk,
629                                cdims);
630                  }                  }
631              }              }
632          }          }
# Line 696  Line 654 
654      Dims[1] = ncb * adims[1];      Dims[1] = ncb * adims[1];
655                                  /* find number of row blocks */                                  /* find number of row blocks */
656      for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];      for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];
657      Dims[0] = nrb * adims[0];      Dims[0] = (nrb + 1) * adims[0]; /* +1 because of 0-based indices */
658      nnz = length(AxP);      nnz = length(AxP);
659    
660      if (nc == 1) {              /* x slot is in the correct order */      if (nc == 1) {              /* x slot is in the correct order */
# Line 747  Line 705 
705      UNPROTECT(1);      UNPROTECT(1);
706      return val;      return val;
707  }  }
708    
709    SEXP dgBCMatrix_to_dgTMatrix(SEXP A)
710    {
711        SEXP val = PROTECT(NEW_OBJECT(MAKE_CLASS("dgTMatrix"))),
712            ApP = GET_SLOT(A, Matrix_pSym),
713            AxP = GET_SLOT(A, Matrix_xSym);
714        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)), *Ap = INTEGER(ApP),
715            *bdims = INTEGER(GET_SLOT(val, Matrix_DimSym)),
716            *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
717            i, j, k, kk, ncb = length(ApP) - 1, nnz = length(AxP), nrb;
718        int *Aj = expand_column_pointers(ncb, Ap, Calloc(nnz, int)),
719            *Bi = INTEGER(ALLOC_SLOT(val, Matrix_iSym, INTSXP, nnz)),
720            *Bj = INTEGER(ALLOC_SLOT(val, Matrix_jSym, INTSXP, nnz)),
721            nblk = adims[2], nc = adims[1], nr = adims[0];
722        int sz = nc * nr;
723        int *ai = Calloc(sz, int), *aj = Calloc(sz, int);
724        double *Ax = REAL(AxP),
725            *Bx = REAL(ALLOC_SLOT(val, Matrix_xSym, REALSXP, nnz));
726    
727        Memcpy(Bx, Ax, nnz);        /* x slot stays as is but w/o dim attribute */
728    
729        bdims[1] = ncb * adims[1];  /* bdims - need to find number of row blocks */
730        for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];
731        bdims[0] = (nrb + 1) * adims[0]; /* +1 because of 0-based indices */
732    
733        for (j = 0; j < nc; j++) {  /* arrays of inner indices */
734            for (i = 0; i < nr; i++) {
735                int ind = j * nc + i;
736                ai[ind] = i;
737                aj[ind] = j;
738            }
739        }
740        for (i = 0, k = 0; k < nblk; k++) {
741            for (kk = 0; kk < sz; kk++) {
742                Bi[k * sz + kk] = Ai[k] * nr + ai[kk];
743                Bj[k * sz + kk] = Aj[k] * nc + aj[kk];
744            }
745        }
746    
747        Free(Aj); Free(ai); Free(aj);
748        UNPROTECT(1);
749        return val;
750    }

Legend:
Removed from v.535  
changed lines
  Added in v.536

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