SCM

SCM Repository

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

Diff of /pkg/src/dgBCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 303, Wed Oct 27 21:33:03 2004 UTC revision 304, Wed Oct 27 21:35:13 2004 UTC
# Line 7  Line 7 
7          xp = GET_SLOT(x, Matrix_xSym),          xp = GET_SLOT(x, Matrix_xSym),
8          dp = getAttrib(xp, R_DimSymbol);          dp = getAttrib(xp, R_DimSymbol);
9      int *pv = INTEGER(pp),      int *pv = INTEGER(pp),
         *iv = INTEGER(ip),  
10          *dim = INTEGER(dp),          *dim = INTEGER(dp),
11          ncol = length(pp) - 1;          ncol = length(pp) - 1;
12      int nnz = pv[ncol];      int nnz = pv[ncol];
13    
14      if (!isReal(xp))      if (!(isReal(xp) && isArray(xp)))
15          return ScalarString(mkChar("slot x should be a real array"));          return ScalarString(mkChar("slot x should be a real array"));
16      if (length(dp) != 3)      if (length(dp) != 3)
17          return ScalarString(mkChar("slot x should be a 3-dimensional array"));          return ScalarString(mkChar("slot x should be a 3-dimensional array"));
# Line 46  Line 45 
45   * @param b matrix to be multiplied   * @param b matrix to be multiplied
46   * @param ldb leading dimension of b as declared in the calling   * @param ldb leading dimension of b as declared in the calling
47   *        routine   *        routine
48     * @param beta scalar multiplier of c
49   * @param c product matrix to be modified   * @param c product matrix to be modified
50   * @param ldc leading dimension of b as declared in the calling   * @param ldc leading dimension of c as declared in the calling
51   *        routine   *        routine
52   */   */
53  void  void
# Line 58  Line 58 
58                const double b[], int ldb,                const double b[], int ldb,
59                double beta, double c[], int ldc)                double beta, double c[], int ldc)
60  {  {
61      int lside = (side == 'L' || side == 'l') ? 1 : 0;      int j, kk, lside = (side == 'L' || side == 'l');
62      int tra = (transa == 'T' || transa == 't') ? 1 : 0;      int ncb, nrb, sz = nr * nc, tra = (transa == 'T' || transa == 't');
63    
64      if (lside && !tra) {        /* only case written so far */      if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0)
65          int j, ncb = k/nc, nrb = m/nr, sz = nc * nr;          error("improper dims m=%d, n=%d, k=%d, nr=%d, nc=%d",
   
         if (nr < 1 || nc < 1 || m < 0 || n < 0 || k < 0 ||  
             m % nr || k % nc)  
             error("incompatible dims m=%d, n=%d, k=%d, nr=%d, nc=%d",  
66                    m, n, k, nr, nc);                    m, n, k, nr, nc);
         if (ldb < k) error("incompatible dims k=%d, ldb=%d", k, ldb);  
67          if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);          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          for (j = 0; j < ncb; j++) {          for (j = 0; j < ncb; j++) {
84              int i, i2 = ap[j+1];              int j2 = ap[j + 1];
85              for (i = ap[j]; i < i2; i++) {              for (kk = ap[j]; kk < j2; kk++) {
86                  int rb = ai[i];                  int ii = ai[kk];
87                  if (rb < 0 || rb >= nrb)                  if (ii < 0 || ii >= nrb)
88                      error("incompatible dims m=%d, nr=%d, rb=%d",                      error("improper row index ii=%d, nrb=%d", ii, nrb);
89                            m, nr, rb);                  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,                  F77_CALL(dgemm)("N", "N", &nr, &n, &nc,
96                                  &alpha, ax + i * sz, &nr,                                      &alpha, ax + j * sz, &nr,
97                                  b + j * nc, &ldb,                                  b + j * nc, &ldb,
98                                  &beta, c + rb * nr, &ldc);                                      &beta, c + ii * nr, &ldc);
99                    }
100              }              }
101          }          }
   
102      } else {      } else {
103          error("Call to cscBlocked_mm must have side == 'L' and transa == 'N'");          error("Call to cscBlocked_mm must have side == 'L'");
104        }
105    }
106    
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  }  }

Legend:
Removed from v.303  
changed lines
  Added in v.304

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