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 414, Mon Jan 3 18:55:11 2005 UTC revision 415, Mon Jan 3 19:05:11 2005 UTC
# Line 1  Line 1 
1  #include "cscBlocked.h"  #include "cscBlocked.h"
2    /* TODO
3     *  - code for trans = 'T' in cscb_syrk
4     *  - code for non-trivial cscb_trmm and cscb_ldl
5     */
6  SEXP cscBlocked_validate(SEXP x)  SEXP cscBlocked_validate(SEXP x)
7  {  {
8      SEXP pp = GET_SLOT(x, Matrix_pSym),      SEXP pp = GET_SLOT(x, Matrix_pSym),
# Line 12  Line 15 
15      int nnz = pv[ncol];      int nnz = pv[ncol];
16    
17      if (!(isReal(xp) && isArray(xp)))      if (!(isReal(xp) && isArray(xp)))
18          return ScalarString(mkChar("slot x should be a real array"));          return mkString("slot x should be a real array");
19      if (length(dp) != 3)      if (length(dp) != 3)
20          return ScalarString(mkChar("slot x should be a 3-dimensional array"));          return mkString("slot x should be a 3-dimensional array");
21      if (length(ip) != nnz)      if (length(ip) != nnz)
22          return ScalarString(mkChar("length of slot i does not matck last element of slot p"));          return mkString("length of slot i does not matck last element of slot p");
23      if (dim[2] != nnz)      if (dim[2] != nnz)
24          return ScalarString(mkChar("third dimension of slot x does not match number of nonzeros"));          return
25                mkString("third dimension of slot x does not match number of nonzeros");
26      return ScalarLogical(1);      return ScalarLogical(1);
27  }  }
28    
# Line 34  Line 38 
38   * @param transa 'T' or 't' for transpose.   * @param transa 'T' or 't' for transpose.
39   * @param m number of rows in C   * @param m number of rows in C
40   * @param n number of columns in C   * @param n number of columns in C
41   * @param k number of columns in B if side = 'L', otherwise   * @param k number of rows in B if side = 'L', otherwise
42   *        number of rows in B.   *        number of columns in B.
43   * @param alpha   * @param alpha
44   * @param nr number of rows per block of A   * @param nr number of rows per block of A
45   * @param nc number of columns per block of A   * @param nc number of columns per block of A
# Line 88  Line 92 
92                      error("improper row index ii=%d, nrb=%d", ii, nrb);                      error("improper row index ii=%d, nrb=%d", ii, nrb);
93                  if (tra) {                  if (tra) {
94                      F77_CALL(dgemm)("T", "N", &nc, &n, &nr,                      F77_CALL(dgemm)("T", "N", &nc, &n, &nr,
95                                      &alpha, ax + j*sz, &nr,                                      &alpha, ax + kk * sz, &nr,
96                                      b + ii * nr, &ldb,                                      b + ii * nr, &ldb,
97                                      &beta, c + j * nc, &ldc);                                      &beta, c + j * nc, &ldc);
98                  } else {                  } else {
99                      F77_CALL(dgemm)("N", "N", &nr, &n, &nc,                      F77_CALL(dgemm)("N", "N", &nr, &n, &nc,
100                                      &alpha, ax + j * sz, &nr,                                      &alpha, ax + kk * sz, &nr,
101                                      b + j*nc, &ldb,                                      b + j*nc, &ldb,
102                                      &beta, c + ii * nr, &ldc);                                      &beta, c + ii * nr, &ldc);
103                  }                  }
# Line 105  Line 109 
109  }  }
110    
111  /**  /**
112     * Perform one of the matrix operations
113     *  C := alpha*op(A)*B + beta*C
114     * or
115     *  C := alpha*B*op(A) + beta*C
116     * where A is a compressed, sparse, blocked matrix and
117     * B and C are dense matrices.
118     *
119     * @param side 'L' or 'l' for left, otherwise right
120     * @param transa 'T' or 't' for transpose.
121     * @param m number of rows in C
122     * @param n number of columns in C
123     * @param k number of rows in B if side = 'L', otherwise
124     *        number of columns in B.
125     * @param alpha
126     * @param A pointer to a cscBlocked object
127     * @param B matrix to be multiplied
128     * @param ldb leading dimension of b as declared in the calling
129     *        routine
130     * @param beta scalar multiplier of c
131     * @param C product matrix to be modified
132     * @param ldc leading dimension of c as declared in the calling
133     *        routine
134     */
135    void
136    cscb_mm(char side, char transa, int m, int n, int k,
137            double alpha, SEXP A,
138            const double B[], int ldb,
139            double beta, double C[], int ldc)
140    {
141        int lside = (side == 'L' || side == 'l'),
142            tra = (transa == 'T' || transa == 't');
143        SEXP AxP = GET_SLOT(A, Matrix_xSym),
144            ApP = GET_SLOT(A, Matrix_pSym);
145        int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
146            *Ap = INTEGER(ApP),
147            *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
148            ancb = length(ApP) - 1, /* number of column blocks */
149            anrb;                   /* number of row blocks */
150        int absz = adims[0] * adims[1]; /* block size */
151        int j;
152        double *Ax = REAL(AxP);
153    
154        if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);
155        if (lside) {
156            /* B is of size k by n */
157            if (ldb < k)
158                error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",
159                      k, ldb, n, adims[0], adims[1]);
160            if (tra) {
161                if (m % adims[1] || k % adims[0])
162                    error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",
163                          m, k, adims[0], adims[1]);
164                if (ancb != m/adims[1])
165                    error("incompatible LT dims m=%d, ancb=%d, adims=[%d,%d,%d]",
166                          m, ancb, adims[0], adims[1], adims[3]);
167                anrb = k/adims[0];
168            } else {
169                if (m % adims[0] || k % adims[1])
170                    error("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d",
171                          m, k, adims[0], adims[1]);
172                if (ancb != k/adims[1])
173                    error("incompatible LN dims k=%d, ancb=%d, adims=[%d,%d,%d]",
174                          k, ancb, adims[0], adims[1], adims[3]);
175                anrb = m/adims[0];
176            }
177            for (j = 0; j < ancb; j++) {
178                int kk, j2 = Ap[j + 1];
179                for (kk = Ap[j]; kk < j2; kk++) {
180                    int ii = Ai[kk];
181                    if (ii < 0 || ii >= anrb)
182                        error("improper row index ii=%d, anrb=%d", ii, anrb);
183                    if (tra) {
184                        F77_CALL(dgemm)("T", "N", adims+1, &n, adims,
185                                        &alpha, Ax + kk * absz, adims,
186                                        B + ii * adims[0], &ldb,
187                                        &beta, C + j * adims[1], &ldc);
188                    } else {
189                        F77_CALL(dgemm)("N", "N", adims, &n, adims+1,
190                                        &alpha, Ax + kk * absz, adims,
191                                        B + j*adims[1], &ldb,
192                                        &beta, C + ii * adims[1], &ldc);
193                    }
194                }
195            }
196        } else {
197            /* B is of size m by k */
198            error("Call to cscb_mm must have side == 'L'");
199        }
200    }
201    
202    /**
203   * Invert a triangular sparse blocked matrix.  This is not done in   * Invert a triangular sparse blocked matrix.  This is not done in
204   * place because the number of non-zero blocks in A-inverse can be   * place because the number of non-zero blocks in A-inverse can be
205   * different than the number of non-zero blocks in A.   * different than the number of non-zero blocks in A.
206   *   *
207   * @param upper 'U' indicates upper triangular, 'L' lower   * @param upper 'U' indicates upper triangular, 'L' lower
208   * @param unit 'U' indicates unit diagonal, 'N' non-unit   * @param unit 'U' indicates unit diagonal, 'N' non-unit
209   * @param n number of columns of blocks in A and A-inverse   * @param A Pointer to a triangular cscBlocked object
210   * @param nr number of rows per block in A and A-inverse   * @param Ai Pointer to a triangular cscBlocked object
  * @param nc number of columns per block in A and A-inverse  
  * @param ap vector of column pointers for A  
  * @param ai vector of row indices for A  
  * @param ax contents of the non-zero blocks of A  
  * @param aip vector of column pointers for A-inverse  
  * @param aii vector of row indices for A-inverse  
  * @param aix contents of the non-zero blocks of A-inverse  
211   */   */
212  void  void
213  cscBlocked_tri(char upper, char unit, int n, int nr, int nc,  cscb_tri(char upper, char unit, SEXP A, const int Parent[], SEXP AI)
                const int ap[], const int ai[], const double ax[],  
                int aip[], int aii[], double aix[])  
214  {  {
215      int /* iup = (upper == 'U' || upper == 'u') ? 1 : 0, */      SEXP ApP = GET_SLOT(A, Matrix_pSym),
216          iunit = (unit == 'U' || unit == 'u') ? 1 : 0;          AxP = GET_SLOT(A, Matrix_xSym),
217            AIpP = GET_SLOT(AI, Matrix_pSym),
218            AIxP = GET_SLOT(AI, Matrix_xSym);
219        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
220            *Ap = INTEGER(ApP),
221            *AIi = INTEGER(GET_SLOT(AI, Matrix_iSym)),
222            *AIp = INTEGER(AIpP),
223            *adim = INTEGER(getAttrib(AxP, R_DimSymbol)),
224            *aidim = INTEGER(getAttrib(AxP, R_DimSymbol)),
225            ancb = length(ApP) - 1,
226            aincb = length(AIpP) - 1,
227            iup = (upper == 'U' || upper == 'u'),
228            iunit = (unit == 'U' || unit == 'u');
229    
230        if (aidim[0] != adim[0] || aidim[1] != adim[1] || aidim[2] != adim[2] ||
231            adim[0] != adim[1] || ancb != aincb)
232            error("Incompatible dimensions A[%d,%d,%d]x%d; AI[%d,%d,%d]x%d",
233                  adim[0], adim[1], adim[2], ancb,
234                  aidim[0], aidim[1], aidim[2], aincb);
235      if (!iunit)      if (!iunit)
236          error("Code for non-unit triangular matrices not yet written");          error("Code for non-unit triangular matrices not yet written");
237      if (ap[n] > 0)      if (adim[2] > 0)
238          error("Code for non-trivial unit inverse not yet written");          error("Code for non-trivial unit inverse not yet written");
     else {  
         if (aip[n] == 0) return;  
         error ("Structure of A and A-inverse does not agree");  
     }  
239  }  }
240    
241  /**  /**
# Line 173  Line 272 
272   * @param beta scalar multiplier of c   * @param beta scalar multiplier of c
273   * @param C compressed sparse blocked symmetric matrix to be updated   * @param C compressed sparse blocked symmetric matrix to be updated
274   */   */
275  SEXP cscBlocked_syrk(SEXP uplo, SEXP trans, SEXP alpha, SEXP A, SEXP beta, SEXP C)  void cscb_syrk(char uplo, char trans, double alpha, SEXP A, double beta, SEXP C)
276  {  {
277      SEXP axp = GET_SLOT(A, Matrix_xSym),      SEXP AxP = GET_SLOT(A, Matrix_xSym),
278          app = GET_SLOT(A, Matrix_pSym),          ApP = GET_SLOT(A, Matrix_pSym),
279          cxp = GET_SLOT(C, Matrix_xSym),          CxP = GET_SLOT(C, Matrix_xSym),
280          cpp = GET_SLOT(C, Matrix_pSym);          CpP = GET_SLOT(C, Matrix_pSym);
281      int *adims = INTEGER(getAttrib(axp, R_DimSymbol)),      int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
282          *ai = INTEGER(GET_SLOT(A, Matrix_iSym)),          *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
283          *ap = INTEGER(app),          *Ap = INTEGER(ApP),
284          *cdims = INTEGER(getAttrib(cxp, R_DimSymbol)),          *cdims = INTEGER(getAttrib(CxP, R_DimSymbol)),
285          *ci = INTEGER(GET_SLOT(A, Matrix_iSym)),          *Ci = INTEGER(GET_SLOT(C, Matrix_iSym)),
286          *cp = INTEGER(cpp),          *Cp = INTEGER(CpP),
287            iup = (uplo == 'U' || uplo == 'u'),
288            itr = (trans == 'T' || trans == 't'),
289          j, k,          j, k,
290          nca = length(app) - 1,          ancb = length(ApP) - 1,
291          ncc = length(cpp) - 1,          cncb = length(CpP) - 1;
292          npairs,      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;
293          inplace;      int scalar = (adims[0] == 1 && adims[1] == 1),
294      char *ul, *tr;          asz = adims[0] * adims[1],
295      double *ax = REAL(axp),          csz = cdims[0] * cdims[1];
296          *cx = REAL(cxp),  
297          bta = asReal(beta);  
   
     if (length(uplo) != 1 || length(trans) != 1 || length(alpha) != 1 ||  
         length(beta) != 1 || !isReal(alpha) || !isReal(beta) ||  
         !isString(uplo) || !isString(trans))  
         error("uplo and trans must be character scalars, alpha and beta real scalars");  
298      if (cdims[0] != cdims[1]) error("blocks in C must be square");      if (cdims[0] != cdims[1]) error("blocks in C must be square");
299      ul = CHAR(STRING_ELT(uplo, 0));      if (itr) {
     tr = CHAR(STRING_ELT(trans, 0));  
     if (toupper(tr[0]) == 'T')  
300          error("Code for trans == 'T' not yet written");          error("Code for trans == 'T' not yet written");
301  /* FIXME: Write the other version */                                  /* FIXME: Write this part */
302      if (adims[0] != cdims[0])      } else {
303          error("Dimension inconsistency in blocks: dim(A)=[%d,,], dim(C)=[%d,,]",          if (adims[1] != cdims[0])
304                adims[0], cdims[0]);              error("Inconsistent dimensions: A[%d,%d,%d]x%d, C[%d,%d,%d]x%d",
305                      adims[0], adims[1], adims[2], ancb,
306                      cdims[0], cdims[1], cdims[2], cncb);
307                                  /* check the row indices */                                  /* check the row indices */
308      for (k = 0; k < adims[2]; k++) {      for (k = 0; k < adims[2]; k++) {
309          int aik = ai[k];              int aik = Ai[k];
310          if (aik < 0 || aik >= ncc)              if (aik < 0 || aik >= cncb)
311              error("Row index %d = %d is out of range [0, %d]",              error("Row index %d = %d is out of range [0, %d]",
312                    k, ai[k], ncc - 1);                        k, aik, cncb - 1);
313            }
314                                    /* multiply C by beta */
315            if (beta != 1.)
316                for (j = 0; j < csz * cdims[2]; j++) Cx[j] *= beta;
317                                    /* individual products */
318            for (j = 0; j < ancb; j++) {
319                int k, kk, k2 = Ap[j+1];
320                for (k = Ap[j]; k < k2; k++) {
321                    int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);
322    
323                    if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
324                    else F77_CALL(dsyrk)(&uplo, "N", cdims, adims + 1,
325                                         &alpha, Ax + k * asz, adims,
326                                         &one, Cx + K * csz, cdims);
327    
328                    for (kk = k+1; kk < k2; kk++) {
329                        int jj = Ai[kk];
330                        K = (iup) ? Ind(jj, ii, Cp, Ci) : Ind(ii, jj, Cp, Ci);
331                        if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];
332                        else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
333                                             &alpha, Ax + ((iup)?kk:k) * asz, adims,
334                                             Ax + ((iup)?k:kk) * asz, adims,
335                                             &one, Cx + K * asz, cdims);
336      }      }
                                 /* check if C can be updated in place */  
     inplace = 1;  
     npairs = 0;  
     for (j = 0; j < nca; j++) {  
         int k, kk, k2 = ap[j+1];  
         int nnz = k2 - ap[j];  
         npairs += (nnz * (nnz - 1)) / 2;  
         for (k = ap[j]; k < k2; k++) {  
 /* FIXME: This check assumes uplo == 'L' */  
             for (kk = k; kk < k2; kk++) {  
                 if (Ind(ai[k], ai[kk], cp, ci) < 0) {  
                     inplace = 0;  
                     break;  
337                  }                  }
338              }              }
             if (!inplace) break;  
339          }          }
340      }      }
                                 /* multiply C by beta */  
     for (j = 0; j < cdims[0]*cdims[1]*cdims[2]; j++) cx[j] *= bta;  
     if (inplace) {  
         int scalar = (adims[0] == 1 && adims[1] == 1);  
341    
342          for (j = 0; j < nca; j++) {  /**
343              int k, kk, k2 = ap[j+1];   * Create the LDL' decomposition of the positive definite symmetric
344              for (k = ap[j]; k < k2; k++) {   * cscBlocked matrix A (upper triangle stored) in L and D.
345                  int ii = ai[k];   *
346                  for (kk = k; kk < k2; kk++) {   * @param A pointer to a cscBlocked object containing the upper
347                      int jj = ai[kk], K = Ind(ii, jj, cp, ci);   * triangle of a positive definite symmetric matrix.
348                      if (scalar) cx[K] += ax[k] * ax[kk];   * @param Parent the parent array for A
349                      else {   * @param L pointer to a cscBlocked object to hold L
350     * @param D pointer to a 3D array to hold D
351     *
352     * @return n the number of column blocks in A for success.  A value
353     * less than n indicates the first column block whose diagonal was not
354     * positive definite.
355     */
356    int
357    cscb_ldl(SEXP A, const int Parent[], SEXP L, SEXP D)
358    {
359        SEXP ApP = GET_SLOT(A, Matrix_pSym),
360            AxP = GET_SLOT(A, Matrix_xSym);
361        int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
362            diag, j, n = length(ApP) - 1;
363        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
364            *Ap = INTEGER(ApP),
365            *Li = INTEGER(GET_SLOT(L, Matrix_iSym)),
366            *Lp = INTEGER(GET_SLOT(L, Matrix_pSym));
367        double *Lx = REAL(GET_SLOT(L, Matrix_xSym)),
368            *Ax = REAL(AxP), *Dx = REAL(D);
369    
370        for (j = 0, diag = 1; j < n; j++) { /* check for trivial structure */
371            if (Parent[j] >= 0) {diag = 0; break;}
372        }
373        if (diag) {
374            Memcpy(REAL(D), Ax, adims[0] * adims[1] * adims[2]);
375            return n;
376        }
377        if (adims[0] == 1 && adims[1] == 1) {
378            return R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent, Li, Lx, Dx,
379                                 (int *) NULL, (int *) NULL);
380        } else {
381    
382            error("code for nontrivial blocked L not yet written");
383                      }                      }
384        return -1;                  /* keep -Wall happy */
385                  }                  }
386    
387    /**
388     * Perform one of the cscBlocked-matrix operations B := alpha*op(A)*B
389     * or B := alpha*B*op(A)
390     *
391     * @param side 'L' or 'R' for left or right
392     * @param uplo 'U' or 'L' for upper or lower
393     * @param transa 'T' or 'N' for transpose or no transpose
394     * @param diag 'U' or 'N' for unit diagonal or non-unit
395     * @param A pointer to a triangular cscBlocked object
396     * @param B pointer to the contents of the matrix B
397     * @param m number of rows in B
398     * @param n number of columns in B
399     * @param ldb leading dimension of B as declared in the calling function
400     */
401    void cscb_trmm(char side, char uplo, char transa, char diag,
402                   double alpha, SEXP A, double B[], int m, int n, int ldb)
403    {
404        int ileft = (side == 'L' || side == 'l'),
405            iup = (uplo == 'U' || uplo == 'u'),
406            itr = (transa == 'T' || transa == 't'),
407            iunit = (diag == 'U' || diag == 'u');
408        SEXP ApP = GET_SLOT(A, Matrix_pSym),
409            AxP = GET_SLOT(A, Matrix_xSym);
410        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
411            *Ap = INTEGER(ApP),
412            *xdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
413            i, j, nb = length(ApP) - 1;
414    
415        if (xdims[0] != xdims[1])
416            error("Argument A to cscb_trmm is not triangular");
417        if (alpha != 1.0) {
418            for (j = 0; j < n; j++) { /* scale by alpha */
419                for (i = 0; i < m; i++)
420                    B[i + j * ldb] *= alpha;
421            }
422        }
423        if (iunit && xdims[2] < 1) return; /* A is the identity */
424        error("Code for non-identity cases of cscb_trmm not yet written");
425        iup += 0;                   /* keep -Wall happy */
426        ileft += 0;
427        itr += 0;
428              }              }
429    
430    /**
431     * Perform one of the operations B := alpha*op(A)*B or
432     * B := alpha*B*op(A) where A and B are both cscBlocked.
433     *
434     * @param side 'L' or 'R' for left or right
435     * @param uplo 'U' or 'L' for upper or lower
436     * @param transa 'T' or 'N' for transpose or no transpose
437     * @param diag 'U' or 'N' for unit diagonal or non-unit
438     * @param alpha scalar multiplier
439     * @param A pointer to a triangular cscBlocked object
440     * @param B pointer to a general cscBlocked matrix
441     */
442    void cscb_trcbm(char side, char uplo, char transa, char diag,
443                    double alpha, SEXP A, SEXP B)
444    {
445        int ileft = (side == 'L' || side == 'l'),
446            iup = (uplo == 'U' || uplo == 'u'),
447            itr = (transa == 'T' || transa == 't'),
448            iunit = (diag == 'U' || diag == 'u');
449        SEXP ApP = GET_SLOT(A, Matrix_pSym),
450            AxP = GET_SLOT(A, Matrix_xSym),
451            BpP = GET_SLOT(B, Matrix_pSym),
452            BxP = GET_SLOT(B, Matrix_xSym);
453        int *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
454            *Ap = INTEGER(ApP),
455            *Bi = INTEGER(GET_SLOT(B, Matrix_iSym)),
456            *Bp = INTEGER(BpP),
457            *axdims = INTEGER(getAttrib(AxP, R_DimSymbol)),
458            *bxdims = INTEGER(getAttrib(BxP, R_DimSymbol)),
459            ncbA = length(ApP) - 1,
460            ncbB = length(BpP) - 1;
461        int i, nbx = bxdims[0] * bxdims[1] * bxdims[2];
462    
463        if (axdims[0] != axdims[1])
464            error("Argument A to cscb_trcbm is not triangular");
465        if (alpha != 1.0) {
466            for (i = 0; i < nbx; i++) { /* scale by alpha */
467                REAL(BxP)[i] *= alpha;
468          }          }
469      }      }
470      return C;      if (iunit && axdims[2] < 1) return; /* A is the identity */
471        error("Code for non-trivial cscb_trcbm not yet written");
472  }  }
473    
474    /**
475     * Perform one of the matrix-matrix operations
476     *      C := alpha*op(A)*op(B) + beta*C
477     * on compressed, sparse, blocked matrices.
478     *
479     * @param transa 'T' for transpose of A, else 'N'
480     * @param transb 'T' for transpose of B, else 'N'
481     * @param alpha scalar multiplier
482     * @param A pointer to a cscBlocked object
483     * @param B pointer to a cscBlocked object
484     * @param beta scalar multiplier
485     * @param C pointer to a cscBlocked object
486     */
487    void cscb_cscbm(char transa, char transb, double alpha, SEXP A, SEXP B,
488                 double beta, SEXP C)
489    {
490        int ta = (transa == 'T' || transa == 't') ? 1 : 0,
491            tb = (transb == 'T' || transb == 't') ? 1 : 0;
492        SEXP ApP = GET_SLOT(A, Matrix_pSym),
493            AxP = GET_SLOT(A, Matrix_xSym),
494            BpP = GET_SLOT(B, Matrix_pSym),
495            BxP = GET_SLOT(B, Matrix_xSym),
496            CxP = GET_SLOT(C, Matrix_xSym);
497        int *Ap = INTEGER(ApP),
498            *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
499            *Bp = INTEGER(BpP),
500            *Bi = INTEGER(GET_SLOT(B, Matrix_iSym)),
501            *Cp = INTEGER(GET_SLOT(C, Matrix_pSym)),
502            *Ci = INTEGER(GET_SLOT(C, Matrix_iSym)),
503            *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
504            *bdims = INTEGER(getAttrib(BxP, R_DimSymbol)),
505            *cdims = INTEGER(getAttrib(CxP, R_DimSymbol)),
506            nca = length(ApP) - 1,
507            ncb = length(BpP) - 1;
508        int ablk = adims[0] * adims[1],
509            bblk = bdims[0] * bdims[1],
510            cblk = cdims[0] * cdims[1];
511        double *Ax = REAL(AxP),
512            *Bx = REAL(BxP),
513            *Cx = REAL(CxP),
514            one = 1.0;
515    
516        if ((!ta) && tb) {          /* transposed crossproduct */
517            int jj;
518    
519            if (adims[1] != bdims[1] ||
520                adims[0] != cdims[0] ||
521                bdims[0] != cdims[1])
522                error("C[%d,%d,%d] := A[%d,%d,%d] %*% t(B[%d,%d,%d])",
523                      cdims[0], cdims[1], cdims[2],
524                      adims[0], adims[1], adims[2],
525                      bdims[0], bdims[1], bdims[2]);
526            if (nca != ncb)
527                error("C := A(ncblocks = %d) %*% t(B(ncblocks = %d)", nca, ncb);
528            if (beta != 1.) {       /* scale C by beta */
529                int ctot = cdims[0] * cdims[1] * cdims[2];
530                for (jj = 0; jj < ctot; jj++) Cx[jj] *= beta;
531            }
532            for (jj = 0; jj < nca; jj++) {
533                int ia, ib, a2 = Ap[jj + 1], b2 = Bp[jj + 1];
534                for (ia = Ap[jj]; ia < a2; ia++) {
535                    for (ib = Bp[jj]; ib < b2; ib++) {
536                        int cind = Ind(Ai[ia], Bi[ib], Cp, Ci);
537                        if (cind < 0)
538                            error("Invalid index [%d,%d]", Ai[ia], Bi[ib]);
539                        F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,
540                                        &alpha, Ax + ia * ablk, adims,
541                                        Bx + ib * bblk, bdims, &one,
542                                        Cx + cind * cblk, cdims);
543                    }
544                }
545            }
546            return;
547        }
548        error("Code not yet written");
549    }

Legend:
Removed from v.414  
changed lines
  Added in v.415

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