SCM

SCM Repository

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

Annotation of /pkg/src/cscBlocked.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 304 - (view) (download) (as text)

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 :     *dim = INTEGER(dp),
11 :     ncol = length(pp) - 1;
12 :     int nnz = pv[ncol];
13 :    
14 : bates 304 if (!(isReal(xp) && isArray(xp)))
15 : bates 298 return ScalarString(mkChar("slot x should be a real array"));
16 :     if (length(dp) != 3)
17 :     return ScalarString(mkChar("slot x should be a 3-dimensional array"));
18 :     if (length(ip) != nnz)
19 :     return ScalarString(mkChar("length of slot i does not matck last element of slot p"));
20 :     if (dim[2] != nnz)
21 :     return ScalarString(mkChar("third dimension of slot x does not match number of nonzeros"));
22 :     return ScalarLogical(1);
23 :     }
24 :    
25 : bates 300 /**
26 :     * Perform one of the matrix operations
27 :     * C := alpha*op(A)*B + beta*C
28 :     * or
29 :     * C := alpha*B*op(A) + beta*C
30 :     * where A is a compressed, sparse, blocked matrix and
31 :     * B and C are dense matrices.
32 :     *
33 :     * @param side 'L' or 'l' for left, otherwise right
34 :     * @param transa 'T' or 't' for transpose.
35 :     * @param m number of rows in C
36 :     * @param n number of columns in C
37 :     * @param k number of columns in B if side = 'L', otherwise
38 :     * number of rows in B.
39 :     * @param alpha
40 :     * @param nr number of rows per block of A
41 :     * @param nc number of columns per block of A
42 :     * @param ap vector of column pointers in A
43 :     * @param ai vector of row indices in A
44 :     * @param ax contents of non-zero blocks of A
45 :     * @param b matrix to be multiplied
46 :     * @param ldb leading dimension of b as declared in the calling
47 :     * routine
48 : bates 304 * @param beta scalar multiplier of c
49 : bates 300 * @param c product matrix to be modified
50 : bates 304 * @param ldc leading dimension of c as declared in the calling
51 : bates 300 * 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 : bates 304 int j, kk, lside = (side == 'L' || side == 'l');
62 :     int ncb, nrb, sz = nr * nc, tra = (transa == 'T' || transa == 't');
63 : bates 300
64 : bates 304 if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)
65 :     error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",
66 : bates 300 m, n, k, nr, nc);
67 : bates 304 if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);
68 :     if (lside) {
69 :     if (ldb < k)
70 :     error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",
71 :     k, ldb, n, nr, nc);
72 :     if (tra) {
73 :     if (m % nc || k % nr)
74 :     error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",
75 :     m, k, nr, nc);
76 :     nrb = k/nr; ncb = m/nc;
77 :     } else {
78 :     if (m % nr || k % nc)
79 :     error("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d",
80 :     m, k, nr, nc);
81 :     nrb = m/nr; ncb = k/nc;
82 :     }
83 : bates 300 for (j = 0; j < ncb; j++) {
84 : bates 304 int j2 = ap[j + 1];
85 :     for (kk = ap[j]; kk < j2; kk++) {
86 :     int ii = ai[kk];
87 :     if (ii < 0 || ii >= nrb)
88 :     error("improper row index ii=%d, nrb=%d", ii, nrb);
89 :     if (tra) {
90 :     F77_CALL(dgemm)("T", "N", &nc, &n, &nr,
91 :     &alpha, ax + j*sz, &nr,
92 :     b + ii * nr, &ldb,
93 :     &beta, c + j * nc, &ldc);
94 :     } else {
95 :     F77_CALL(dgemm)("N", "N", &nr, &n, &nc,
96 :     &alpha, ax + j * sz, &nr,
97 :     b + j*nc, &ldb,
98 :     &beta, c + ii * nr, &ldc);
99 :     }
100 : bates 300 }
101 :     }
102 :     } else {
103 : bates 304 error("Call to cscBlocked_mm must have side == 'L'");
104 : bates 300 }
105 :     }
106 : bates 304
107 :     /**
108 :     * Invert a triangular sparse blocked matrix. This is not done in
109 :     * place because the number of non-zero blocks in A-inverse can be
110 :     * different than the number of non-zero blocks in A.
111 :     *
112 :     * @param upper 'U' indicates upper triangular, 'L' lower
113 :     * @param unit 'U' indicates unit diagonal, 'N' non-unit
114 :     * @param n number of columns of blocks in A and A-inverse
115 :     * @param nr number of rows per block in A and A-inverse
116 :     * @param nc number of columns per block in A and A-inverse
117 :     * @param ap vector of column pointers for A
118 :     * @param ai vector of row indices for A
119 :     * @param ax contents of the non-zero blocks of A
120 :     * @param aip vector of column pointers for A-inverse
121 :     * @param aii vector of row indices for A-inverse
122 :     * @param aix contents of the non-zero blocks of A-inverse
123 :     */
124 :     void
125 :     cscBlocked_tri(char upper, char unit, int n, int nr, int nc,
126 :     const int ap[], const int ai[], const double ax[],
127 :     int aip[], int aii[], double aix[])
128 :     {
129 :     int iup = (upper == 'U' || upper == 'u') ? 1 : 0,
130 :     iunit = (unit == 'U' || unit == 'u') ? 1 : 0;
131 :    
132 :     if (!iunit)
133 :     error("Code for non-unit triangular matrices not yet written");
134 :     if (ap[n] > 0)
135 :     error("Code for non-trivial unit inverse not yet written");
136 :     else {
137 :     if (aip[n] == 0) return;
138 :     error ("Structure of A and A-inverse does not agree");
139 :     }
140 :     }

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