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 415 - (view) (download) (as text)
Original Path: pkg/src/cscBlocked.c

1 : bates 298 #include "cscBlocked.h"
2 : bates 415 /* TODO
3 :     * - code for trans = 'T' in cscb_syrk
4 :     * - code for non-trivial cscb_trmm and cscb_ldl
5 :     */
6 : bates 298 SEXP cscBlocked_validate(SEXP x)
7 :     {
8 :     SEXP pp = GET_SLOT(x, Matrix_pSym),
9 :     ip = GET_SLOT(x, Matrix_iSym),
10 :     xp = GET_SLOT(x, Matrix_xSym),
11 :     dp = getAttrib(xp, R_DimSymbol);
12 :     int *pv = INTEGER(pp),
13 :     *dim = INTEGER(dp),
14 :     ncol = length(pp) - 1;
15 :     int nnz = pv[ncol];
16 :    
17 : bates 304 if (!(isReal(xp) && isArray(xp)))
18 : bates 415 return mkString("slot x should be a real array");
19 : bates 298 if (length(dp) != 3)
20 : bates 415 return mkString("slot x should be a 3-dimensional array");
21 : bates 298 if (length(ip) != nnz)
22 : bates 415 return mkString("length of slot i does not matck last element of slot p");
23 : bates 298 if (dim[2] != nnz)
24 : bates 415 return
25 :     mkString("third dimension of slot x does not match number of nonzeros");
26 : bates 298 return ScalarLogical(1);
27 :     }
28 :    
29 : bates 300 /**
30 :     * Perform one of the matrix operations
31 :     * C := alpha*op(A)*B + beta*C
32 :     * or
33 :     * C := alpha*B*op(A) + beta*C
34 :     * where A is a compressed, sparse, blocked matrix and
35 :     * B and C are dense matrices.
36 :     *
37 :     * @param side 'L' or 'l' for left, otherwise right
38 :     * @param transa 'T' or 't' for transpose.
39 :     * @param m number of rows in C
40 :     * @param n number of columns in C
41 : bates 415 * @param k number of rows in B if side = 'L', otherwise
42 :     * number of columns in B.
43 : bates 300 * @param alpha
44 :     * @param nr number of rows per block of A
45 :     * @param nc number of columns per block of A
46 :     * @param ap vector of column pointers in A
47 :     * @param ai vector of row indices in A
48 :     * @param ax contents of non-zero blocks of A
49 :     * @param b matrix to be multiplied
50 :     * @param ldb leading dimension of b as declared in the calling
51 :     * routine
52 : bates 304 * @param beta scalar multiplier of c
53 : bates 300 * @param c product matrix to be modified
54 : bates 304 * @param ldc leading dimension of c as declared in the calling
55 : bates 300 * routine
56 :     */
57 :     void
58 :     cscBlocked_mm(char side, char transa, int m, int n, int k,
59 :     double alpha, int nr, int nc,
60 :     const int ap[], const int ai[],
61 :     const double ax[],
62 :     const double b[], int ldb,
63 :     double beta, double c[], int ldc)
64 :     {
65 : bates 304 int j, kk, lside = (side == 'L' || side == 'l');
66 :     int ncb, nrb, sz = nr * nc, tra = (transa == 'T' || transa == 't');
67 : bates 300
68 : bates 304 if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)
69 :     error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",
70 : bates 300 m, n, k, nr, nc);
71 : bates 304 if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);
72 :     if (lside) {
73 :     if (ldb < k)
74 :     error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",
75 :     k, ldb, n, nr, nc);
76 :     if (tra) {
77 :     if (m % nc || k % nr)
78 :     error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",
79 :     m, k, nr, nc);
80 :     nrb = k/nr; ncb = m/nc;
81 :     } else {
82 :     if (m % nr || k % nc)
83 :     error("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d",
84 :     m, k, nr, nc);
85 :     nrb = m/nr; ncb = k/nc;
86 :     }
87 : bates 300 for (j = 0; j < ncb; j++) {
88 : bates 304 int j2 = ap[j + 1];
89 :     for (kk = ap[j]; kk < j2; kk++) {
90 :     int ii = ai[kk];
91 :     if (ii < 0 || ii >= nrb)
92 :     error("improper row index ii=%d, nrb=%d", ii, nrb);
93 :     if (tra) {
94 :     F77_CALL(dgemm)("T", "N", &nc, &n, &nr,
95 : bates 415 &alpha, ax + kk * sz, &nr,
96 : bates 304 b + ii * nr, &ldb,
97 :     &beta, c + j * nc, &ldc);
98 :     } else {
99 :     F77_CALL(dgemm)("N", "N", &nr, &n, &nc,
100 : bates 415 &alpha, ax + kk * sz, &nr,
101 : bates 304 b + j*nc, &ldb,
102 :     &beta, c + ii * nr, &ldc);
103 :     }
104 : bates 300 }
105 :     }
106 :     } else {
107 : bates 304 error("Call to cscBlocked_mm must have side == 'L'");
108 : bates 300 }
109 :     }
110 : bates 304
111 :     /**
112 : bates 415 * 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 : bates 304 * 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
205 :     * different than the number of non-zero blocks in A.
206 :     *
207 :     * @param upper 'U' indicates upper triangular, 'L' lower
208 :     * @param unit 'U' indicates unit diagonal, 'N' non-unit
209 : bates 415 * @param A Pointer to a triangular cscBlocked object
210 :     * @param Ai Pointer to a triangular cscBlocked object
211 : bates 304 */
212 :     void
213 : bates 415 cscb_tri(char upper, char unit, SEXP A, const int Parent[], SEXP AI)
214 : bates 304 {
215 : bates 415 SEXP ApP = GET_SLOT(A, Matrix_pSym),
216 :     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 : bates 304 if (!iunit)
236 :     error("Code for non-unit triangular matrices not yet written");
237 : bates 415 if (adim[2] > 0)
238 : bates 304 error("Code for non-trivial unit inverse not yet written");
239 :     }
240 : bates 339
241 :     /**
242 :     * Search for the element in a compressed sparse matrix at a given row and column
243 :     *
244 :     * @param row row index
245 :     * @param col column index
246 :     * @param cp column pointers
247 :     * @param ci row indices
248 :     *
249 :     * @return index of element in ci, if it exists, else -1
250 :     */
251 :     static R_INLINE
252 :     int Ind(int row, int col, const int cp[], const int ci[])
253 :     {
254 :     int i, i2 = cp[col + 1];
255 :     for (i = cp[col]; i < i2; i++)
256 :     if (ci[i] == row) return i;
257 :     return -1;
258 :     }
259 :    
260 :     /**
261 :     * Perform one of the matrix operations
262 :     * C := alpha*A*A' + beta*C,
263 :     * or
264 :     * C := alpha*A'*A + beta*C,
265 :     * where A is a compressed, sparse, blocked matrix and
266 :     * C is a compressed, sparse, symmetric blocked matrix.
267 :     *
268 :     * @param uplo 'U' or 'u' for upper triangular storage, else lower.
269 :     * @param trans 'T' or 't' for transpose.
270 :     * @param alpha scalar multiplier of outer product
271 :     * @param A compressed sparse blocked matrix
272 :     * @param beta scalar multiplier of c
273 :     * @param C compressed sparse blocked symmetric matrix to be updated
274 :     */
275 : bates 415 void cscb_syrk(char uplo, char trans, double alpha, SEXP A, double beta, SEXP C)
276 : bates 339 {
277 : bates 415 SEXP AxP = GET_SLOT(A, Matrix_xSym),
278 :     ApP = GET_SLOT(A, Matrix_pSym),
279 :     CxP = GET_SLOT(C, Matrix_xSym),
280 :     CpP = GET_SLOT(C, Matrix_pSym);
281 :     int *adims = INTEGER(getAttrib(AxP, R_DimSymbol)),
282 :     *Ai = INTEGER(GET_SLOT(A, Matrix_iSym)),
283 :     *Ap = INTEGER(ApP),
284 :     *cdims = INTEGER(getAttrib(CxP, R_DimSymbol)),
285 :     *Ci = INTEGER(GET_SLOT(C, Matrix_iSym)),
286 :     *Cp = INTEGER(CpP),
287 :     iup = (uplo == 'U' || uplo == 'u'),
288 :     itr = (trans == 'T' || trans == 't'),
289 : bates 339 j, k,
290 : bates 415 ancb = length(ApP) - 1,
291 :     cncb = length(CpP) - 1;
292 :     double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;
293 :     int scalar = (adims[0] == 1 && adims[1] == 1),
294 :     asz = adims[0] * adims[1],
295 :     csz = cdims[0] * cdims[1];
296 : bates 339
297 : bates 415
298 : bates 339 if (cdims[0] != cdims[1]) error("blocks in C must be square");
299 : bates 415 if (itr) {
300 : bates 339 error("Code for trans == 'T' not yet written");
301 : bates 415 /* FIXME: Write this part */
302 :     } else {
303 :     if (adims[1] != cdims[0])
304 :     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 : bates 339 /* check the row indices */
308 : bates 415 for (k = 0; k < adims[2]; k++) {
309 :     int aik = Ai[k];
310 :     if (aik < 0 || aik >= cncb)
311 :     error("Row index %d = %d is out of range [0, %d]",
312 :     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 : bates 339 }
337 :     }
338 :     }
339 :     }
340 : bates 415 }
341 :    
342 :     /**
343 :     * Create the LDL' decomposition of the positive definite symmetric
344 :     * cscBlocked matrix A (upper triangle stored) in L and D.
345 :     *
346 :     * @param A pointer to a cscBlocked object containing the upper
347 :     * triangle of a positive definite symmetric matrix.
348 :     * @param Parent the parent array for A
349 :     * @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 :     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 : bates 339 }
544 :     }
545 :     }
546 : bates 415 return;
547 : bates 339 }
548 : bates 415 error("Code not yet written");
549 : bates 339 }

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