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 421, Mon Jan 3 20:10:42 2005 UTC revision 422, Tue Jan 4 23:45:15 2005 UTC
# Line 151  Line 151 
151      int j;      int j;
152      double *Ax = REAL(AxP);      double *Ax = REAL(AxP);
153    
154      if (ldc < n) error("incompatible dims n=%d, ldc=%d", n, ldc);      if (ldc < m) error("incompatible dims m=%d, ldc=%d", m, ldc);
155      if (lside) {      if (lside) {
156          /* B is of size k by n */          /* B is of size k by n */
157          if (ldb < k)          if (ldb < k)
# Line 189  Line 189 
189                      F77_CALL(dgemm)("N", "N", adims, &n, adims+1,                      F77_CALL(dgemm)("N", "N", adims, &n, adims+1,
190                                      &alpha, Ax + kk * absz, adims,                                      &alpha, Ax + kk * absz, adims,
191                                      B + j*adims[1], &ldb,                                      B + j*adims[1], &ldb,
192                                      &beta, C + ii * adims[1], &ldc);                                      &beta, C + ii * adims[0], &ldc);
193                  }                  }
194              }              }
195          }          }
# Line 286  Line 286 
286          *Cp = INTEGER(CpP),          *Cp = INTEGER(CpP),
287          iup = (uplo == 'U' || uplo == 'u'),          iup = (uplo == 'U' || uplo == 'u'),
288          itr = (trans == 'T' || trans == 't'),          itr = (trans == 'T' || trans == 't'),
289          j, k,          j, k;
         ancb = length(ApP) - 1,  
         cncb = length(CpP) - 1;  
290      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;      double *Ax = REAL(AxP), *Cx = REAL(CxP), one = 1.;
291      int scalar = (adims[0] == 1 && adims[1] == 1),      int scalar = (adims[0] == 1 && adims[1] == 1),
292          asz = adims[0] * adims[1],          asz = adims[0] * adims[1],
# Line 300  Line 298 
298          error("Code for trans == 'T' not yet written");          error("Code for trans == 'T' not yet written");
299                                  /* FIXME: Write this part */                                  /* FIXME: Write this part */
300      } else {      } else {
301          if (adims[1] != cdims[0])          if (adims[0] != cdims[0])
302              error("Inconsistent dimensions: A[%d,%d,%d]x%d, C[%d,%d,%d]x%d",              error("Inconsistent dimensions: A[%d,%d,%d], C[%d,%d,%d]",
303                    adims[0], adims[1], adims[2], ancb,                    adims[0], adims[1], adims[2],
304                    cdims[0], cdims[1], cdims[2], cncb);                    cdims[0], cdims[1], cdims[2]);
305                                  /* check the row indices */                                  /* check the row indices */
306          for (k = 0; k < adims[2]; k++) {          for (k = 0; k < adims[2]; k++) {
307              int aik = Ai[k];              int aik = Ai[k];
308              if (aik < 0 || aik >= cncb)              if (aik < 0 || aik >= cdims[2])
309                  error("Row index %d = %d is out of range [0, %d]",                  error("Row index %d = %d is out of range [0, %d]",
310                        k, aik, cncb - 1);                        k, aik, cdims[2] - 1);
311          }          }
312                                  /* multiply C by beta */                                  /* multiply C by beta */
313          if (beta != 1.)          if (beta != 1.)
314              for (j = 0; j < csz * cdims[2]; j++) Cx[j] *= beta;              for (j = 0; j < csz * cdims[2]; j++) Cx[j] *= beta;
315                                  /* individual products */                                  /* individual products */
316          for (j = 0; j < ancb; j++) {          for (j = 0; j < adims[2]; j++) {
317              int k, kk, k2 = Ap[j+1];              int k, kk, k2 = Ap[j+1];
318              for (k = Ap[j]; k < k2; k++) {              for (k = Ap[j]; k < k2; k++) {
319                  int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);                  int ii = Ai[k], K = Ind(ii, ii, Cp, Ci);

Legend:
Removed from v.421  
changed lines
  Added in v.422

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