SCM

SCM Repository

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

Annotation of /pkg/src/dgBCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 300 - (view) (download) (as text)
Original Path: pkg/src/cscBlocked.c

1 : bates 298 #include "cscBlocked.h"
2 :    
3 :     SEXP cscBlocked_validate(SEXP x)
4 :     {
5 :     SEXP pp = GET_SLOT(x, Matrix_pSym),
6 :     ip = GET_SLOT(x, Matrix_iSym),
7 :     xp = GET_SLOT(x, Matrix_xSym),
8 :     dp = getAttrib(xp, R_DimSymbol);
9 :     int *pv = INTEGER(pp),
10 :     *iv = INTEGER(ip),
11 :     *dim = INTEGER(dp),
12 :     ncol = length(pp) - 1;
13 :     int nnz = pv[ncol];
14 :    
15 :     if (!isReal(xp))
16 :     return ScalarString(mkChar("slot x should be a real array"));
17 :     if (length(dp) != 3)
18 :     return ScalarString(mkChar("slot x should be a 3-dimensional array"));
19 :     if (length(ip) != nnz)
20 :     return ScalarString(mkChar("length of slot i does not matck last element of slot p"));
21 :     if (dim[2] != nnz)
22 :     return ScalarString(mkChar("third dimension of slot x does not match number of nonzeros"));
23 :     return ScalarLogical(1);
24 :     }
25 :    
26 : bates 300 /**
27 :     * Perform one of the matrix operations
28 :     * C := alpha*op(A)*B + beta*C
29 :     * or
30 :     * C := alpha*B*op(A) + beta*C
31 :     * where A is a compressed, sparse, blocked matrix and
32 :     * B and C are dense matrices.
33 :     *
34 :     * @param side 'L' or 'l' for left, otherwise right
35 :     * @param transa 'T' or 't' for transpose.
36 :     * @param m number of rows in C
37 :     * @param n number of columns in C
38 :     * @param k number of columns in B if side = 'L', otherwise
39 :     * number of rows in B.
40 :     * @param alpha
41 :     * @param nr number of rows per block of A
42 :     * @param nc number of columns per block of A
43 :     * @param ap vector of column pointers in A
44 :     * @param ai vector of row indices in A
45 :     * @param ax contents of non-zero blocks of A
46 :     * @param b matrix to be multiplied
47 :     * @param ldb leading dimension of b as declared in the calling
48 :     * routine
49 :     * @param c product matrix to be modified
50 :     * @param ldc leading dimension of b as declared in the calling
51 :     * routine
52 :     */
53 :     void
54 :     cscBlocked_mm(char side, char transa, int m, int n, int k,
55 :     double alpha, int nr, int nc,
56 :     const int ap[], const int ai[],
57 :     const double ax[],
58 :     const double b[], int ldb,
59 :     double beta, double c[], int ldc)
60 :     {
61 :     int lside = (side == 'L' || side == 'l') ? 1 : 0;
62 :     int tra = (transa == 'T' || transa == 't') ? 1 : 0;
63 :    
64 :     if (lside && !tra) { /* only case written so far */
65 :     int j, ncb = k/nc, nrb = m/nr, sz = nc * nr;
66 :    
67 :     if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0 ||
68 :     m % nr || k % nc)
69 :     error("incompatible dims m=%d, n=%d, k=%d, nr=%d, nc=%d",
70 :     m, n, k, nr, nc);
71 :     if (ldb < k) error("incompatible dims k=%d, ldb=%d", k, ldb);
72 :     if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);
73 :     for (j = 0; j < ncb; j++) {
74 :     int i, i2 = ap[j+1];
75 :     for (i = ap[j]; i < i2; i++) {
76 :     int rb = ai[i];
77 :     if (rb < 0 || rb >= nrb)
78 :     error("incompatible dims m=%d, nr=%d, rb=%d",
79 :     m, nr, rb);
80 :     F77_CALL(dgemm)("N", "N", &nr, &n, &nc,
81 :     &alpha, ax + i * sz, &nr,
82 :     b + j * nc, &ldb,
83 :     &beta, c + rb * nr, &ldc);
84 :     }
85 :     }
86 :    
87 :     } else {
88 :     error("Call to cscBlocked_mm must have side == 'L' and transa == 'N'");
89 :     }
90 :     }

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