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 446, Wed Jan 19 19:08:01 2005 UTC revision 447, Fri Jan 21 12:17:56 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
5     *  - replace calls to Ind by calls to Tind (unless the error checking is different)
6   */   */
7    
8  SEXP cscBlocked_validate(SEXP x)  SEXP cscBlocked_validate(SEXP x)
9  {  {
10      SEXP pp = GET_SLOT(x, Matrix_pSym),      SEXP pp = GET_SLOT(x, Matrix_pSym),
# Line 77  Line 79 
79   *        routine   *        routine
80   */   */
81  void  void
82  cscBlocked_mm(char side, char transa, int m, int n, int k,  cscBlocked_mm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa, int m, int n, int k,
83                double alpha, int nr, int nc,                double alpha, int nr, int nc,
84                const int ap[], const int ai[],                const int ap[], const int ai[],
85                const double ax[],                const double ax[],
86                const double b[], int ldb,                const double b[], int ldb,
87                double beta, double c[], int ldc)                double beta, double c[], int ldc)
88  {  {
89      int j, kk, lside = (side == 'L' || side == 'l');      int j, kk;
90      int ncb, nrb, sz = nr * nc, tra = (transa == 'T' || transa == 't');      int ncb, nrb, sz = nr * nc;
91    
92      if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)      if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)
93          error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",          error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",
94                    m, n, k, nr, nc);                    m, n, k, nr, nc);
95      if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);      if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);
96      if (lside) {      if (side == LFT) {
97          if (ldb < k)          if (ldb < k)
98              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",
99                    k, ldb, n, nr, nc);                    k, ldb, n, nr, nc);
100          if (tra) {          if (transa == TRN) {
101              if (m % nc || k % nr)              if (m % nc || k % nr)
102                  error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",                  error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",
103                        m, k, nr, nc);                        m, k, nr, nc);
# Line 112  Line 114 
114                  int ii = ai[kk];                  int ii = ai[kk];
115                  if (ii < 0 || ii >= nrb)                  if (ii < 0 || ii >= nrb)
116                      error("improper row index ii=%d, nrb=%d", ii, nrb);                      error("improper row index ii=%d, nrb=%d", ii, nrb);
117                  if (tra) {                  if (transa == TRN) {
118                      F77_CALL(dgemm)("T", "N", &nc, &n, &nr,                      F77_CALL(dgemm)("T", "N", &nc, &n, &nr,
119                                      &alpha, ax + kk * sz, &nr,                                      &alpha, ax + kk * sz, &nr,
120                                      b + ii * nr, &ldb,                                      b + ii * nr, &ldb,
# Line 155  Line 157 
157   *        routine   *        routine
158   */   */
159  void  void
160  cscb_mm(char side, char transa, int m, int n, int k,  cscb_mm(enum CBLAS_SIDE side, enum CBLAS_TRANSPOSE transa,
161          double alpha, SEXP A,          int m, int n, int k, double alpha, SEXP A,
162          const double B[], int ldb,          const double B[], int ldb, double beta, double C[], int ldc)
         double beta, double C[], int ldc)  
163  {  {
164      int lside = (side == 'L' || side == 'l'),      int lside = (side == LFT);
         tra = (transa == 'T' || transa == 't');  
165      SEXP AxP = GET_SLOT(A, Matrix_xSym),      SEXP AxP = GET_SLOT(A, Matrix_xSym),
166          ApP = GET_SLOT(A, Matrix_pSym);          ApP = GET_SLOT(A, Matrix_pSym);
167      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
# Line 179  Line 179 
179          if (ldb < k)          if (ldb < k)
180              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",
181                    k, ldb, n, adims[0], adims[1]);                    k, ldb, n, adims[0], adims[1]);
182          if (tra) {          if (transa == TRN) {
183              if (m % adims[1] || k % adims[0])              if (m % adims[1] || k % adims[0])
184                  error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",                  error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",
185                        m, k, adims[0], adims[1]);                        m, k, adims[0], adims[1]);
186              if (ancb != m/adims[1])              if (ancb != m/adims[1])
187                  error("incompatible LT dims m=%d, ancb=%d, adims=[%d,%d,%d]",                  error("incompatible LT dims m=%d, ancb=%d, adims=[%d,%d,%d]",
188                        m, ancb, adims[0], adims[1], adims[3]);                        m, ancb, adims[0], adims[1], adims[2]);
189              anrb = k/adims[0];              anrb = k/adims[0];
190          } else {          } else {
191              if (m % adims[0] || k % adims[1])              if (m % adims[0] || k % adims[1])
# Line 193  Line 193 
193                        m, k, adims[0], adims[1]);                        m, k, adims[0], adims[1]);
194              if (ancb != k/adims[1])              if (ancb != k/adims[1])
195                  error("incompatible LN dims k=%d, ancb=%d, adims=[%d,%d,%d]",                  error("incompatible LN dims k=%d, ancb=%d, adims=[%d,%d,%d]",
196                        k, ancb, adims[0], adims[1], adims[3]);                        k, ancb, adims[0], adims[1], adims[2]);
197              anrb = m/adims[0];              anrb = m/adims[0];
198          }          }
199          for (j = 0; j < ancb; j++) {          for (j = 0; j < ancb; j++) {
# Line 202  Line 202 
202                  int ii = Ai[kk];                  int ii = Ai[kk];
203                  if (ii < 0 || ii >= anrb)                  if (ii < 0 || ii >= anrb)
204                      error("improper row index ii=%d, anrb=%d", ii, anrb);                      error("improper row index ii=%d, anrb=%d", ii, anrb);
205                  if (tra) {                  if (transa == TRN) {
206                      F77_CALL(dgemm)("T", "N", adims+1, &n, adims,                      F77_CALL(dgemm)("T", "N", adims+1, &n, adims,
207                                      &alpha, Ax + kk * absz, adims,                                      &alpha, Ax + kk * absz, adims,
208                                      B + ii * adims[0], &ldb,                                      B + ii * adims[0], &ldb,
# Line 222  Line 222 
222  }  }
223    
224  /**  /**
  * Invert a triangular sparse blocked matrix.  This is not done in  
  * place because the number of non-zero blocks in A-inverse can be  
  * different than the number of non-zero blocks in A.  
  *  
  * @param upper 'U' indicates upper triangular, 'L' lower  
  * @param unit 'U' indicates unit diagonal, 'N' non-unit  
  * @param A Pointer to a triangular cscBlocked object  
  * @param Ai Pointer to a triangular cscBlocked object  
  */  
 void  
 cscb_tri(char upper, char unit, SEXP A, const int Parent[], SEXP AI)  
 {  
     SEXP ApP = GET_SLOT(A, Matrix_pSym),  
         AxP = GET_SLOT(A, Matrix_xSym),  
         AIpP = GET_SLOT(AI, Matrix_pSym),  
         AIxP = GET_SLOT(AI, Matrix_xSym);  
     int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),  
         *Ap = INTEGER(ApP),  
         *AIi = INTEGER(GET_SLOT(AI, Matrix_iSym)),  
         *AIp = INTEGER(AIpP),  
         *adim = INTEGER(getAttrib(AxP, R_DimSymbol)),  
         *aidim = INTEGER(getAttrib(AxP, R_DimSymbol)),  
         ancb = length(ApP) - 1,  
         aincb = length(AIpP) - 1,  
         iup = (upper == 'U' || upper == 'u'),  
         iunit = (unit == 'U' || unit == 'u');  
   
     if (aidim[0] != adim[0] || aidim[1] != adim[1] || aidim[2] != adim[2] ||  
         adim[0] != adim[1] || ancb != aincb)  
         error("Incompatible dimensions A[%d,%d,%d]x%d; AI[%d,%d,%d]x%d",  
               adim[0], adim[1], adim[2], ancb,  
               aidim[0], aidim[1], aidim[2], aincb);  
     if (!iunit)  
         error("Code for non-unit triangular matrices not yet written");  
     if (adim[2] > 0)  
         error("Code for non-trivial unit inverse not yet written");  
 }  
   
 /**  
225   * Search for the element in a compressed sparse matrix at a given row and column   * Search for the element in a compressed sparse matrix at a given row and column
226   *   *
227   * @param row row index   * @param row row index
# Line 294  Line 255 
255   * @param beta scalar multiplier of c   * @param beta scalar multiplier of c
256   * @param C compressed sparse blocked symmetric matrix to be updated   * @param C compressed sparse blocked symmetric matrix to be updated
257   */   */
258  void cscb_syrk(char uplo, char trans, double alpha, SEXP A, double beta, SEXP C)  void
259    cscb_syrk(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE trans,
260              double alpha, SEXP A,
261              double beta, SEXP C)
262  {  {
263      SEXP AxP = GET_SLOT(A, Matrix_xSym),      SEXP AxP = GET_SLOT(A, Matrix_xSym),
264          ApP = GET_SLOT(A, Matrix_pSym),          ApP = GET_SLOT(A, Matrix_pSym),
# Line 343  Line 307 
307    
308                  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);
309                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
310                  else F77_CALL(dsyrk)(&uplo, "N", cdims, adims + 1,                  else F77_CALL(dsyrk)((uplo == UPP)?"U":"L", "N", cdims, adims + 1,
311                                       &alpha, Ax + k * asz, adims,                                       &alpha, Ax + k * asz, adims,
312                                       &one, Cx + K * csz, cdims);                                       &one, Cx + K * csz, cdims);
   
313                  for (kk = k+1; kk < k2; kk++) {                  for (kk = k+1; kk < k2; kk++) {
314                      int jj = Ai[kk];                      int jj = Ai[kk];
315                      K = (iup) ? Ind(ii, jj, Cp, Ci) : Ind(jj, ii, Cp, Ci);                      K = (iup) ? Ind(ii, jj, Cp, Ci) : Ind(jj, ii, Cp, Ci);
# Line 493  Line 456 
456   * @param n number of columns in B   * @param n number of columns in B
457   * @param ldb leading dimension of B as declared in the calling function   * @param ldb leading dimension of B as declared in the calling function
458   */   */
459  void cscb_trmm(char side, char uplo, char transa, char diag,  void
460    cscb_trmm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo,
461              enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
462                 double alpha, SEXP A, double B[], int m, int n, int ldb)                 double alpha, SEXP A, double B[], int m, int n, int ldb)
463  {  {
464      int ileft = (side == 'L' || side == 'l'),      int ileft = (side == LFT),
465          iup = (uplo == 'U' || uplo == 'u'),          iup = (uplo == UPP),
466          itr = (transa == 'T' || transa == 't'),          itr = (transa == TRN),
467          iunit = (diag == 'U' || diag == 'u');          iunit = (diag == UNT);
468      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
469          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
470      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
# Line 535  Line 500 
500   * @param n number of columns in B   * @param n number of columns in B
501   * @param ldb leading dimension of B as declared in the calling function   * @param ldb leading dimension of B as declared in the calling function
502   */   */
503  void cscb_trsm(char uplo, char transa, char diag,  void
504    cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
505                 double alpha, SEXP A, double B[], int m, int n, int ldb)                 double alpha, SEXP A, double B[], int m, int n, int ldb)
506  {  {
507      int iup = (uplo == 'U' || uplo == 'u'),      int iup = (uplo == UPP),
508          itr = (transa == 'T' || transa == 't'),          itr = (transa == TRN),
509          iunit = (diag == 'U' || diag == 'u');          iunit = (diag == UNT);
510      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
511          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
512      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),      int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
# Line 593  Line 559 
559   * @param A pointer to a triangular cscBlocked object   * @param A pointer to a triangular cscBlocked object
560   * @param B pointer to a general cscBlocked matrix   * @param B pointer to a general cscBlocked matrix
561   */   */
562  void cscb_trcbm(char side, char uplo, char transa, char diag,  void
563    cscb_trcbm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
564                  double alpha, SEXP A, SEXP B)                  double alpha, SEXP A, SEXP B)
565  {  {
566      int ileft = (side == 'L' || side == 'l'),      int ileft = (side == LFT),
567          iup = (uplo == 'U' || uplo == 'u'),          iup = (uplo == UPP),
568          itr = (transa == 'T' || transa == 't'),          itr = (transa == TRN),
569          iunit = (diag == 'U' || diag == 'u');          iunit = (diag == UNT);
570      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
571          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
572          BpP = GET_SLOT(B, Matrix_pSym),          BpP = GET_SLOT(B, Matrix_pSym),
# Line 660  Line 627 
627   * @param A pointer to a triangular cscBlocked object   * @param A pointer to a triangular cscBlocked object
628   * @param B pointer to a general cscBlocked matrix   * @param B pointer to a general cscBlocked matrix
629   */   */
630  void cscb_trcbsm(char side, char uplo, char transa, char diag,  void
631    cscb_trcbsm(enum CBLAS_SIDE side, enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,
632                   double alpha, SEXP A, const int Parent[], SEXP B)                   double alpha, SEXP A, const int Parent[], SEXP B)
633  {  {
634      int ileft = (side == 'L' || side == 'l'),      int ileft = (side == LFT),
635          iup = (uplo == 'U' || uplo == 'u'),          iup = (uplo == UPP),
636          itr = (transa == 'T' || transa == 't'),          itr = (transa == TRN),
637          iunit = (diag == 'U' || diag == 'u');          iunit = (diag == UNT);
638      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
639          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
640          BpP = GET_SLOT(B, Matrix_pSym),          BpP = GET_SLOT(B, Matrix_pSym),
# Line 734  Line 702 
702   * @param beta scalar multiplier   * @param beta scalar multiplier
703   * @param C pointer to a cscBlocked object   * @param C pointer to a cscBlocked object
704   */   */
705  void cscb_cscbm(char transa, char transb, double alpha, SEXP A, SEXP B,  void
706               double beta, SEXP C)  cscb_cscbm(enum CBLAS_TRANSPOSE transa, enum CBLAS_TRANSPOSE transb,
707               double alpha, SEXP A, SEXP B, double beta, SEXP C)
708  {  {
     int ta = (transa == 'T' || transa == 't') ? 1 : 0,  
         tb = (transb == 'T' || transb == 't') ? 1 : 0;  
709      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
710          AxP = GET_SLOT(A, Matrix_xSym),          AxP = GET_SLOT(A, Matrix_xSym),
711          BpP = GET_SLOT(B, Matrix_pSym),          BpP = GET_SLOT(B, Matrix_pSym),
# Line 763  Line 730 
730          *Cx = REAL(CxP),          *Cx = REAL(CxP),
731          one = 1.0;          one = 1.0;
732    
733      if ((!ta) && tb) {          /* transposed crossproduct */      if ((transa == NTR) && transb == TRN) {             /* transposed crossproduct */
734          int jj;          int jj;
735    
736          if (adims[1] != bdims[1] ||          if (adims[1] != bdims[1] ||

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

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