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 299, Fri Oct 8 14:48:06 2004 UTC revision 300, Tue Oct 12 15:44:15 2004 UTC
# Line 23  Line 23 
23      return ScalarLogical(1);      return ScalarLogical(1);
24  }  }
25    
26    /**
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    }

Legend:
Removed from v.299  
changed lines
  Added in v.300

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