SCM

SCM Repository

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

Diff of /pkg/src/cscBlocked.c

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

revision 447, Fri Jan 21 12:17:56 2005 UTC revision 448, Fri Jan 21 21:59:45 2005 UTC
# Line 2  Line 2 
2  /* TODO  /* TODO
3   *  - code for trans = 'T' in cscb_syrk   *  - code for trans = 'T' in cscb_syrk
4   *  - code for non-trivial cscb_trmm and cscb_ldl   *  - code for non-trivial cscb_trmm and cscb_ldl
  *  - replace calls to Ind by calls to Tind (unless the error checking is different)  
5   */   */
6    
7  SEXP cscBlocked_validate(SEXP x)  SEXP cscBlocked_validate(SEXP x)
# Line 65  Line 64 
64   * @param k number of rows in B if side = 'L', otherwise   * @param k number of rows in B if side = 'L', otherwise
65   *        number of columns in B.   *        number of columns in B.
66   * @param alpha   * @param alpha
  * @param nr number of rows per block of A  
  * @param nc number of columns per block of A  
  * @param ap vector of column pointers in A  
  * @param ai vector of row indices in A  
  * @param ax contents of non-zero blocks of A  
  * @param b matrix to be multiplied  
  * @param ldb leading dimension of b as declared in the calling  
  *        routine  
  * @param beta scalar multiplier of c  
  * @param c product matrix to be modified  
  * @param ldc leading dimension of c as declared in the calling  
  *        routine  
  */  
 void  
 cscBlocked_mm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m, int n, int k,  
               double alpha, int nr, int nc,  
               const int ap[], const int ai[],  
               const double ax[],  
               const double b[], int ldb,  
               double beta, double c[], int ldc)  
 {  
     int j, kk;  
     int ncb, nrb, sz = nr * nc;  
   
     if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)  
         error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",  
                   m, n, k, nr, nc);  
     if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);  
     if (side == LFT) {  
         if (ldb < k)  
             error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",  
                   k, ldb, n, nr, nc);  
         if (transa == TRN) {  
             if (m % nc || k % nr)  
                 error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",  
                       m, k, nr, nc);  
             nrb = k/nr; ncb = m/nc;  
         } else {  
             if (m % nr || k % nc)  
                 error("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d",  
                       m, k, nr, nc);  
             nrb = m/nr; ncb = k/nc;  
         }  
         for (j = 0; j < ncb; j++) {  
             int j2 = ap[j + 1];  
             for (kk = ap[j]; kk < j2; kk++) {  
                 int ii = ai[kk];  
                 if (ii < 0 || ii >= nrb)  
                     error("improper row index ii=%d, nrb=%d", ii, nrb);  
                 if (transa == TRN) {  
                     F77_CALL(dgemm)("T", "N", &nc, &n, &nr,  
                                     &alpha, ax + kk * sz, &nr,  
                                     b + ii * nr, &ldb,  
                                     &beta, c + j * nc, &ldc);  
                 } else {  
                     F77_CALL(dgemm)("N", "N", &nr, &n, &nc,  
                                     &alpha, ax + kk * sz, &nr,  
                                     b + j*nc, &ldb,  
                                     &beta, c + ii * nr, &ldc);  
                 }  
             }  
         }  
     } else {  
         error("Call to cscBlocked_mm must have side == 'L'");  
     }  
 }  
   
 /**  
  * Perform one of the matrix operations  
  *  C := alpha*op(A)*B + beta*C  
  * or  
  *  C := alpha*B*op(A) + beta*C  
  * where A is a compressed, sparse, blocked matrix and  
  * B and C are dense matrices.  
  *  
  * @param side 'L' or 'l' for left, otherwise right  
  * @param transa 'T' or 't' for transpose.  
  * @param m number of rows in C  
  * @param n number of columns in C  
  * @param k number of rows in B if side = 'L', otherwise  
  *        number of columns in B.  
  * @param alpha  
67   * @param A pointer to a cscBlocked object   * @param A pointer to a cscBlocked object
68   * @param B matrix to be multiplied   * @param B matrix to be multiplied
69   * @param ldb leading dimension of b as declared in the calling   * @param ldb leading dimension of b as declared in the calling
# Line 161  Line 78 
78          int m, int n, int k, double alpha, SEXP A,          int m, int n, int k, double alpha, SEXP A,
79          const double B[], int ldb, double beta, double C[], int ldc)          const double B[], int ldb, double beta, double C[], int ldc)
80  {  {
     int lside = (side == LFT);  
81      SEXP AxP = GET_SLOT(A, Matrix_xSym),      SEXP AxP = GET_SLOT(A, Matrix_xSym),
82          ApP = GET_SLOT(A, Matrix_pSym);          ApP = GET_SLOT(A, Matrix_pSym);
83      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
# Line 174  Line 90 
90      double *Ax = REAL(AxP);      double *Ax = REAL(AxP);
91    
92      if (ldc < m) error("incompatible dims m=%d, ldc=%d", m, ldc);      if (ldc < m) error("incompatible dims m=%d, ldc=%d", m, ldc);
93      if (lside) {      if (side == LFT) {
94          /* B is of size k by n */          /* B is of size k by n */
95          if (ldb < k)          if (ldb < k)
96              error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",              error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",
# Line 222  Line 138 
138  }  }
139    
140  /**  /**
  * Search for the element in a compressed sparse matrix at a given row and column  
  *  
  * @param row row index  
  * @param col column index  
  * @param cp column pointers  
  * @param ci row indices  
  *  
  * @return index of element in ci, if it exists, else -1  
  */  
 static R_INLINE  
 int Ind(int row, int col, const int cp[], const int ci[])  
 {  
     int i, i2 = cp[col + 1];  
     for (i = cp[col]; i < i2; i++)  
         if (ci[i] == row) return i;  
     return -1;  
 }  
   
 /**  
141   * Perform one of the matrix operations   * Perform one of the matrix operations
142   *  C := alpha*A*A' + beta*C,   *  C := alpha*A*A' + beta*C,
143   * or   * or
# Line 270  Line 167 
167          *cdims = INTEGER(getAttrib(CxP, R_DimSymbol)),          *cdims = INTEGER(getAttrib(CxP, R_DimSymbol)),
168          *Ci = INTEGER(GET_SLOT(C, Matrix_iSym)),          *Ci = INTEGER(GET_SLOT(C, Matrix_iSym)),
169          *Cp = INTEGER(CpP),          *Cp = INTEGER(CpP),
         iup = (uplo == 'U' || uplo == 'u'),  
         itr = (trans == 'T' || trans == 't'),  
170          j, k;          j, k;
171      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;
172      int scalar = (adims[0] == 1 && adims[1] == 1),      int scalar = (adims[0] == 1 && adims[1] == 1),
# Line 281  Line 176 
176    
177    
178      if (cdims[0] != cdims[1]) error("blocks in C must be square");      if (cdims[0] != cdims[1]) error("blocks in C must be square");
179      if (itr) {      if (trans == TRN) {
180          error("Code for trans == 'T' not yet written");          error("Code for trans == 'T' not yet written");
181                                  /* FIXME: Write this part */                                  /* FIXME: Write this part */
182      } else {      } else {
# Line 303  Line 198 
198          for (j = 0; j < anc; j++) {          for (j = 0; j < anc; j++) {
199              int k, kk, k2 = Ap[j+1];              int k, kk, k2 = Ap[j+1];
200              for (k = Ap[j]; k < k2; k++) {              for (k = Ap[j]; k < k2; k++) {
201                  int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);                  int ii = Ai[k], K = Tind(Ci, Cp, ii, ii);
202    
203                  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);
204                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
# Line 312  Line 207 
207                                       &one, Cx + K * csz, cdims);                                       &one, Cx + K * csz, cdims);
208                  for (kk = k+1; kk < k2; kk++) {                  for (kk = k+1; kk < k2; kk++) {
209                      int jj = Ai[kk];                      int jj = Ai[kk];
210                      K = (iup) ? Ind(ii, jj, Cp, Ci) : Ind(jj, ii, Cp, Ci);                      K = (uplo == UPP) ? Tind(Ci, Cp, ii, jj) : Tind(Ci, Cp, jj, ii);
211    
                     if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, jj);  
212                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];
213                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
214                                           &alpha, Ax + ((iup)?kk:k) * asz, adims,                                           &alpha, Ax + ((uplo==UPP)?kk:k) * asz, adims,
215                                           Ax + ((iup)?k:kk) * asz, adims,                                           Ax + ((uplo==UPP)?k:kk) * asz, adims,
216                                           &one, Cx + K * asz, cdims);                                           &one, Cx + K * asz, cdims);
217                  }                  }
218              }              }
# Line 446  Line 340 
340   * Perform one of the cscBlocked-matrix operations B := alpha*op(A)*B   * Perform one of the cscBlocked-matrix operations B := alpha*op(A)*B
341   * or B := alpha*B*op(A)   * or B := alpha*B*op(A)
342   *   *
343   * @param side 'L' or 'R' for left or right   * @param side
344   * @param uplo 'U' or 'L' for upper or lower   * @param uplo
345   * @param transa 'T' or 'N' for transpose or no transpose   * @param transa
346   * @param diag 'U' or 'N' for unit diagonal or non-unit   * @param diag
347   * @param A pointer to a triangular cscBlocked object   * @param A pointer to a triangular cscBlocked object
348   * @param B pointer to the contents of the matrix B   * @param B contents of the matrix B
349   * @param m number of rows in B   * @param m number of rows in B
350   * @param n number of columns in B   * @param n number of columns in B
351   * @param ldb leading dimension of B as declared in the calling function   * @param ldb leading dimension of B as declared in the calling function
# Line 461  Line 355 
355            enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,            enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
356            double alpha, SEXP A, double B[], int m, int n, int ldb)            double alpha, SEXP A, double B[], int m, int n, int ldb)
357  {  {
     int ileft = (side == LFT),  
         iup = (uplo == UPP),  
         itr = (transa == TRN),  
         iunit = (diag == UNT);  
358      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
359          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
360      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
# Line 480  Line 370 
370                  B[i + j * ldb] *= alpha;                  B[i + j * ldb] *= alpha;
371          }          }
372      }      }
373      if (iunit && xdims[2] < 1) return; /* A is the identity */      if (diag == UNT && xdims[2] < 1) return; /* A is the identity */
374      error("Code for non-identity cases of cscb_trmm not yet written");      error("Code for non-identity cases of cscb_trmm not yet written");
     iup += 0;                   /* keep -Wall happy */  
     ileft += 0;  
     itr += 0;  
375  }  }
376    
377  /**  /**
# Line 504  Line 391 
391  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
392            double alpha, SEXP A, double B[], int m, int n, int ldb)            double alpha, SEXP A, double B[], int m, int n, int ldb)
393  {  {
     int iup = (uplo == UPP),  
         itr = (transa == TRN),  
         iunit = (diag == UNT);  
394      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
395          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
396      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
397          *Ap = INTEGER(ApP),          *Ap = INTEGER(ApP),
398          *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),          *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
399          i, j, nb = length(ApP) - 1;          i, j, nb = length(ApP) - 1;
400      double *Ax = REAL(GET_SLOT(A, Matrix_xSym));      double *Ax = REAL(GET_SLOT(A, Matrix_xSym)), minus1 = -1., one = 1.;
401    
402      if (xdims[0] != xdims[1])      if (xdims[0] != xdims[1])
403          error("Argument A to cscb_trsm is not triangular");          error("Argument A to cscb_trsm is not triangular");
# Line 529  Line 413 
413                  B[i + j * ldb] *= alpha;                  B[i + j * ldb] *= alpha;
414          }          }
415      }      }
416      if (iunit) {      if (diag == UNT) {
417          if (xdims[2] < 1) return; /* A is the identity */          if (xdims[2] < 1) return; /* A is the identity */
418          if (xdims[0] == 1) {          if (xdims[0] == 1) {    /* scalar case */
419              if (iup) error("Code for upper triangle not yet written");              if (uplo == UPP) error("Code for upper triangle not yet written");
420              if (itr) {              if (transa == TRN) {
421                  for (j = 0; j < n; j++)                  for (j = 0; j < n; j++)
422                      R_ldl_ltsolve(m, B + j * ldb, Ap, Ai, Ax);                      R_ldl_ltsolve(m, B + j * ldb, Ap, Ai, Ax);
423              } else {              } else {
# Line 541  Line 425 
425                      R_ldl_lsolve(m, B + j * ldb, Ap, Ai, Ax);                      R_ldl_lsolve(m, B + j * ldb, Ap, Ai, Ax);
426              }              }
427              return;              return;
428            } else {
429                int p, p2, sza = xdims[0] * xdims[0], szb = xdims[0] * n;
430                double *tmp = Calloc(szb, double);
431                if (uplo == UPP) error("Code for upper triangle not yet written");
432                if (transa == TRN) {
433                    for (j = nb - 1; j >= 0; j--) {
434                        p2 = Ap[j+1];
435    
436                        F77_CALL(dlacpy)("A", xdims, &n, B + j * xdims[0], &ldb,
437                                         tmp, xdims);
438                        for (p = Ap[j]; p < p2; p++)
439                            F77_CALL(dgemm)("T", "N", xdims, &n, xdims,
440                                            &minus1, Ax + p * sza, xdims,
441                                            B + Ai[p] * xdims[0], &ldb,
442                                            &one, tmp, xdims);
443                        F77_CALL(dlacpy)("A", xdims, &n, tmp, xdims,
444                                         B + j * xdims[0], &ldb);
445                    }
446                } else {
447                    for (j = 0; j < nb; j++) {
448                        p2 = Ap[j+1];
449    
450                        F77_CALL(dlacpy)("A", xdims, &n, B + j * xdims[0], &ldb,
451                                         tmp, xdims);
452                        for (p = Ap[j]; p < p2; p++)
453                            F77_CALL(dgemm)("N", "N", xdims, &n, xdims,
454                                            &minus1, Ax + p * sza, xdims,
455                                            B + Ai[p] * xdims[0], &ldb,
456                                            &one, tmp, xdims);
457                        F77_CALL(dlacpy)("A", xdims, &n, tmp, xdims,
458                                         B + j * xdims[0], &ldb);
459                    }
460          }          }
         error("Code for non-scalar cscBlocked objects not yet written");  
461      }      }
462      error("Code for non-unit cases of cscb_trsm not yet written");      } else {error("Code for non-unit cases of cscb_trsm not yet written");}
463  }  }
464    
465  /**  /**
466   * Perform one of the operations B := alpha*op(A)*B or   * Perform one of the operations B := alpha*op(A)*B or
467   * B := alpha*B*op(A) where A and B are both cscBlocked.   * B := alpha*B*op(A) where A and B are both cscBlocked.
468   *   *
469   * @param side 'L' or 'R' for left or right   * @param side
470   * @param uplo 'U' or 'L' for upper or lower   * @param uplo
471   * @param transa 'T' or 'N' for transpose or no transpose   * @param transa
472   * @param diag 'U' or 'N' for unit diagonal or non-unit   * @param diag
473   * @param alpha scalar multiplier   * @param alpha scalar multiplier
474   * @param A pointer to a triangular cscBlocked object   * @param A pointer to a triangular cscBlocked object
475   * @param B pointer to a general cscBlocked matrix   * @param B pointer to a general cscBlocked matrix
476   */   */
477  void  void
478  cscb_trcbm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,  cscb_trcbm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo,
479               enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
480             double alpha, SEXP A, SEXP B)             double alpha, SEXP A, SEXP B)
481  {  {
482      int ileft = (side == LFT),      int
         iup = (uplo == UPP),  
         itr = (transa == TRN),  
483          iunit = (diag == UNT);          iunit = (diag == UNT);
484      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
485          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
# Line 588  Line 502 
502              REAL(BxP)[i] *= alpha;              REAL(BxP)[i] *= alpha;
503          }          }
504      }      }
505      if (iunit && axdims[2] < 1) return; /* A is the identity */      if (diag == UNT && axdims[2] < 1) return; /* A is the identity */
506      error("Code for non-trivial cscb_trcbm not yet written");      error("Code for non-trivial cscb_trcbm not yet written");
507  }  }
508    
# Line 628  Line 542 
542   * @param B pointer to a general cscBlocked matrix   * @param B pointer to a general cscBlocked matrix
543   */   */
544  void  void
545  cscb_trcbsm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,  cscb_trcbsm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo,
546                enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
547              double alpha, SEXP A, const int Parent[], SEXP B)              double alpha, SEXP A, const int Parent[], SEXP B)
548  {  {
     int ileft = (side == LFT),  
         iup = (uplo == UPP),  
         itr = (transa == TRN),  
         iunit = (diag == UNT);  
549      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
550          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
551          BpP = GET_SLOT(B, Matrix_pSym),          BpP = GET_SLOT(B, Matrix_pSym),
# Line 657  Line 568 
568              REAL(BxP)[i] *= alpha;              REAL(BxP)[i] *= alpha;
569          }          }
570      }      }
571      if (iunit && axdims[2] < 1) return; /* A is the identity */      if (diag == UNT && axdims[2] < 1) return;   /* A is the identity */
572      if (iunit && axdims[0] == 1) { /* can use R_ldl code */      if (diag == UNT && axdims[0] == 1) { /* can use R_ldl code */
573          if (!ileft && itr) {    /* case required for lmer */          if ((side != LFT) && transa == TRN) {   /* case required for lmer */
574              int *BTp, nnz = bxdims[2], nrbB;              int *BTp, nnz = bxdims[2], nrbB;
575              int *tmp = expand_column_pointers(ncbB, Bp, Calloc(nnz, int));              int *tmp = expand_column_pointers(ncbB, Bp, Calloc(nnz, int));
576              int *BTi = Calloc(nnz, int);              int *BTi = Calloc(nnz, int);
# Line 750  Line 661 
661              int ia, ib, a2 = Ap[jj + 1], b2 = Bp[jj + 1];              int ia, ib, a2 = Ap[jj + 1], b2 = Bp[jj + 1];
662              for (ia = Ap[jj]; ia < a2; ia++) {              for (ia = Ap[jj]; ia < a2; ia++) {
663                  for (ib = Bp[jj]; ib < b2; ib++) {                  for (ib = Bp[jj]; ib < b2; ib++) {
                     int cind = Ind(Ai[ia], Bi[ib], Cp, Ci);  
                     if (cind < 0)  
                         error("Invalid index [%d,%d]", Ai[ia], Bi[ib]);  
664                      F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,                      F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
665                                      &alpha, Ax + ia * ablk, adims,                                      &alpha, Ax + ia * ablk, adims,
666                                      Bx + ib * bblk, bdims, &one,                                      Bx + ib * bblk, bdims, &one,
667                                      Cx + cind * cblk, cdims);                                      Cx + Tind(Ci, Cp, Ai[ia], Bi[ib])*cblk, cdims);
668                  }                  }
669              }              }
670          }          }

Legend:
Removed from v.447  
changed lines
  Added in v.448

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