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 565, Tue Feb 22 01:15:45 2005 UTC revision 582, Mon Feb 28 18:15:21 2005 UTC
# Line 16  Line 16 
16      int nnz = pv[ncol];      int nnz = pv[ncol];
17    
18      if (!(isReal(xp) && isArray(xp)))      if (!(isReal(xp) && isArray(xp)))
19          return mkString("slot x should be a real array");          return mkString(_("slot x should be a real array"));
20      if (length(dp) != 3)      if (length(dp) != 3)
21          return mkString("slot x should be a 3-dimensional array");          return mkString(_("slot x should be a 3-dimensional array"));
22      if (length(ip) != nnz)      if (length(ip) != nnz)
23          return mkString("length of slot i does not matck last element of slot p");          return mkString(_("length of slot i does not match last element of slot p"));
24      if (dim[2] != nnz)      if (dim[2] != nnz)
25          return          return
26              mkString("third dimension of slot x does not match number of nonzeros");              mkString(_("third dimension of slot x does not match number of nonzeros"));
27      return ScalarLogical(1);      return ScalarLogical(1);
28  }  }
29    
# Line 67  Line 67 
67      int j;      int j;
68      double *Ax = REAL(AxP);      double *Ax = REAL(AxP);
69    
70      if (ldc < m) error("incompatible dims m=%d, ldc=%d", m, ldc);      if (ldc < m) error(_("incompatible dims m=%d, ldc=%d"), m, ldc);
71      if (side == LFT) {      if (side == LFT) {
72          /* B is of size k by n */          /* B is of size k by n */
73          if (ldb < k)          if (ldb < k)
74              error("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d",              error(_("incompatible L dims k=%d, ldb=%d, n=%d, nr=%d, nc=%d"),
75                    k, ldb, n, adims[0], adims[1]);                    k, ldb, n, adims[0], adims[1]);
76          if (transa == TRN) {          if (transa == TRN) {
77              if (m % adims[1] || k % adims[0])              if (m % adims[1] || k % adims[0])
78                  error("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d",                  error(_("incompatible LT dims m=%d, k = %d, nr=%d, nc=%d"),
79                        m, k, adims[0], adims[1]);                        m, k, adims[0], adims[1]);
80              if (ancb != m/adims[1])              if (ancb != m/adims[1])
81                  error("incompatible LT dims m=%d, ancb=%d, adims=[%d,%d,%d]",                  error(_("incompatible LT dims m=%d, ancb=%d, adims=[%d,%d,%d]"),
82                        m, ancb, adims[0], adims[1], adims[2]);                        m, ancb, adims[0], adims[1], adims[2]);
83              anrb = k/adims[0];              anrb = k/adims[0];
84          } else {          } else {
85              if (m % adims[0] || k % adims[1])              if (m % adims[0] || k % adims[1])
86                  error("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d",                  error(_("incompatible LN dims m=%d, k = %d, nr=%d, nc=%d"),
87                        m, k, adims[0], adims[1]);                        m, k, adims[0], adims[1]);
88              if (ancb != k/adims[1])              if (ancb != k/adims[1])
89                  error("incompatible LN dims k=%d, ancb=%d, adims=[%d,%d,%d]",                  error(_("incompatible LN dims k=%d, ancb=%d, adims=[%d,%d,%d]"),
90                        k, ancb, adims[0], adims[1], adims[2]);                        k, ancb, adims[0], adims[1], adims[2]);
91              anrb = m/adims[0];              anrb = m/adims[0];
92          }          }
# Line 95  Line 95 
95              for (kk = Ap[j]; kk < j2; kk++) {              for (kk = Ap[j]; kk < j2; kk++) {
96                  int ii = Ai[kk];                  int ii = Ai[kk];
97                  if (ii < 0 || ii >= anrb)                  if (ii < 0 || ii >= anrb)
98                      error("improper row index ii=%d, anrb=%d", ii, anrb);                      error(_("improper row index ii=%d, anrb=%d"), ii, anrb);
99                  if (transa == TRN) {                  if (transa == TRN) {
100                      F77_CALL(dgemm)("T", "N", adims+1, &n, adims,                      F77_CALL(dgemm)("T", "N", adims+1, &n, adims,
101                                      &alpha, Ax + kk * absz, adims,                                      &alpha, Ax + kk * absz, adims,
# Line 111  Line 111 
111          }          }
112      } else {      } else {
113          /* B is of size m by k */          /* B is of size m by k */
114          error("Call to cscb_mm must have side == LFT");          error(_("Call to cscb_mm must have side == LFT"));
115      }      }
116  }  }
117    
# Line 153  Line 153 
153          csz = cdims[0] * cdims[1];          csz = cdims[0] * cdims[1];
154    
155    
156      if (cdims[0] != cdims[1]) error("blocks in C must be square");      if (cdims[0] != cdims[1]) error(_("blocks in C must be square"));
157      if (trans == TRN) {      if (trans == TRN) {
158                                  /* FIXME: Write this part */                                  /* FIXME: Write this part */
159          error("Code for trans == TRN not yet written");          error(_("Code for trans == TRN not yet written"));
160      } else {      } else {
161          if (adims[0] != cdims[0])          if (adims[0] != cdims[0])
162              error("Inconsistent dimensions: A[%d,%d,%d], C[%d,%d,%d]",              error(_("Inconsistent dimensions: A[%d,%d,%d], C[%d,%d,%d]"),
163                    adims[0], adims[1], adims[2],                    adims[0], adims[1], adims[2],
164                    cdims[0], cdims[1], cdims[2]);                    cdims[0], cdims[1], cdims[2]);
165                                  /* check the row indices */                                  /* check the row indices */
166          for (k = 0; k < adims[2]; k++) {          for (k = 0; k < adims[2]; k++) {
167              int aik = Ai[k];              int aik = Ai[k];
168              if (aik < 0 || aik >= cdims[2])              if (aik < 0 || aik >= cdims[2])
169                  error("Row index %d = %d is out of range [0, %d]",                  error(_("Row index %d = %d is out of range [0, %d]"),
170                        k, aik, cdims[2] - 1);                        k, aik, cdims[2] - 1);
171          }          }
172                                  /* multiply C by beta */                                  /* multiply C by beta */
# Line 178  Line 178 
178              for (k = Ap[j]; k < k2; k++) {              for (k = Ap[j]; k < k2; k++) {
179                  int ii = Ai[k], K = check_csc_index(Cp, Ci, ii, ii, 0);                  int ii = Ai[k], K = check_csc_index(Cp, Ci, ii, ii, 0);
180    
181                  if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, ii);                  if (K < 0) error(_("cscb_syrk: C[%d,%d] not defined"), ii, ii);
182                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
183                  else F77_CALL(dsyrk)((uplo == UPP) ? "U" : "L", "N",                  else F77_CALL(dsyrk)((uplo == UPP) ? "U" : "L", "N",
184                                       cdims, adims + 1,                                       cdims, adims + 1,
# Line 244  Line 244 
244          *Ax = REAL(AxP), *Dx = REAL(D), minus1 = -1., one = 1.;          *Ax = REAL(AxP), *Dx = REAL(D), minus1 = -1., one = 1.;
245    
246      if (adims[1] != nci || nci < 1)      if (adims[1] != nci || nci < 1)
247          error("cscb_ldl: dim(A) is [%d, %d, %d]", adims[0], adims[1], adims[2]);          error(_("cscb_ldl: dim(A) is [%d, %d, %d]"), adims[0], adims[1], adims[2]);
248      for (j = 0, diag = 1; j < n; j++) { /* check for trivial structure */      for (j = 0, diag = 1; j < n; j++) { /* check for trivial structure */
249          if (Parent[j] >= 0) {diag = 0; break;}          if (Parent[j] >= 0) {diag = 0; break;}
250      }      }
# Line 279  Line 279 
279              p2 = Ap[k+1];              p2 = Ap[k+1];
280              for (p = Ap[k]; p < p2; p++) {              for (p = Ap[k]; p < p2; p++) {
281                  i = Ai[p];      /* get A[i,k] */                  i = Ai[p];      /* get A[i,k] */
282                  if (i > k) error("cscb_ldl: A has nonzeros below diagonal");                  if (i > k) error(_("cscb_ldl: A has nonzeros below diagonal"));
283                                  /* copy A(i,k) into Y */                                  /* copy A(i,k) into Y */
284                  Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);                  Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);
285                  /* follow path from i to root of etree,                  /* follow path from i to root of etree,
# Line 359  Line 359 
359          i, j/* , nb = length(ApP) - 1 */;          i, j/* , nb = length(ApP) - 1 */;
360    
361      if (xdims[0] != xdims[1])      if (xdims[0] != xdims[1])
362          error("Argument A to cscb_trmm is not triangular");          error(_("Argument A to cscb_trmm is not triangular"));
363      if (alpha != 1.0) {      if (alpha != 1.0) {
364          for (j = 0; j < n; j++) { /* scale by alpha */          for (j = 0; j < n; j++) { /* scale by alpha */
365              for (i = 0; i < m; i++)              for (i = 0; i < m; i++)
# Line 367  Line 367 
367          }          }
368      }      }
369      if (diag == UNT && xdims[2] < 1) return; /* A is the identity */      if (diag == UNT && xdims[2] < 1) return; /* A is the identity */
370      error("Code for non-identity cases of cscb_trmm not yet written");      error(_("Code for non-identity cases of cscb_trmm not yet written"));
371  }  }
372    
373  /**  /**
# Line 398  Line 398 
398      double *Ax = REAL(GET_SLOT(A, Matrix_xSym)), minus1 = -1., one = 1.;      double *Ax = REAL(GET_SLOT(A, Matrix_xSym)), minus1 = -1., one = 1.;
399    
400      if (xdims[0] != xdims[1])      if (xdims[0] != xdims[1])
401          error("Argument A to cscb_trsm is not triangular");          error(_("Argument A to cscb_trsm is not triangular"));
402      if (ldb < m || ldb <= 0 || n <= 0)      if (ldb < m || ldb <= 0 || n <= 0)
403          error("cscb_trsm: inconsistent dims m = %d, n = %d, ldb = %d",          error(_("cscb_trsm: inconsistent dims m = %d, n = %d, ldb = %d"),
404                m, n, ldb);                m, n, ldb);
405      if (m != (nb * xdims[0]))      if (m != (nb * xdims[0]))
406          error("cscb_trsm: inconsistent dims m = %d, A[%d,%d,]x%d",          error(_("cscb_trsm: inconsistent dims m = %d, A[%d,%d,]x%d"),
407                m, xdims[0], xdims[1], xdims[2]);                m, xdims[0], xdims[1], xdims[2]);
408      if (alpha != 1.0) {      if (alpha != 1.0) {
409          for (j = 0; j < n; j++) { /* scale by alpha */          for (j = 0; j < n; j++) { /* scale by alpha */
# Line 414  Line 414 
414      if (diag == UNT) {      if (diag == UNT) {
415          if (xdims[2] < 1) return; /* A is the identity */          if (xdims[2] < 1) return; /* A is the identity */
416          if (xdims[0] == 1) {    /* scalar case */          if (xdims[0] == 1) {    /* scalar case */
417              if (uplo == UPP) error("Code for upper triangle not yet written");              if (uplo == UPP) error(_("Code for upper triangle not yet written"));
418              if (transa == TRN) {              if (transa == TRN) {
419                  for (j = 0; j < n; j++)                  for (j = 0; j < n; j++)
420                      R_ldl_ltsolve(m, B + j * ldb, Ap, Ai, Ax);                      R_ldl_ltsolve(m, B + j * ldb, Ap, Ai, Ax);
# Line 426  Line 426 
426          } else {          } else {
427              int p, p2, sza = xdims[0] * xdims[0];              int p, p2, sza = xdims[0] * xdims[0];
428    
429              if (uplo == UPP) error("Code for upper triangle not yet written");              if (uplo == UPP) error(_("Code for upper triangle not yet written"));
430              if (transa == TRN) {              if (transa == TRN) {
431                  for (j = nb - 1; j >= 0; j--) {                  for (j = nb - 1; j >= 0; j--) {
432                      p2 = Ap[j+1];                      p2 = Ap[j+1];
# Line 447  Line 447 
447                  }                  }
448              }              }
449          }          }
450      } else {error("Code for non-unit cases of cscb_trsm not yet written");}      } else {error(_("Code for non-unit cases of cscb_trsm not yet written"));}
451  }  }
452    
453  /**  /**
# Line 485  Line 485 
485      int i, nbx = bxdims[0] * bxdims[1] * bxdims[2];      int i, nbx = bxdims[0] * bxdims[1] * bxdims[2];
486    
487      if (axdims[0] != axdims[1])      if (axdims[0] != axdims[1])
488          error("Argument A to cscb_trcbm is not triangular");          error(_("Argument A to cscb_trcbm is not triangular"));
489      if (alpha != 1.0) {      if (alpha != 1.0) {
490          for (i = 0; i < nbx; i++) { /* scale by alpha */          for (i = 0; i < nbx; i++) { /* scale by alpha */
491              REAL(BxP)[i] *= alpha;              REAL(BxP)[i] *= alpha;
492          }          }
493      }      }
494      if (diag == UNT && axdims[2] < 1) return; /* A is the identity */      if (diag == UNT && axdims[2] < 1) return; /* A is the identity */
495      error("Code for non-trivial cscb_trcbm not yet written");      error(_("Code for non-trivial cscb_trcbm not yet written"));
496  }  }
497    
498  /**  /**
# Line 529  Line 529 
529      double *Ax = REAL(AxP), *Bx = REAL(BxP);      double *Ax = REAL(AxP), *Bx = REAL(BxP);
530    
531      if (axdims[0] != axdims[1])      if (axdims[0] != axdims[1])
532          error("Argument A to cscb_trcbm is not triangular");          error(_("Argument A to cscb_trcbm is not triangular"));
533      if (alpha != 1.0) {      if (alpha != 1.0) {
534          for (i = 0; i < nbx; i++) { /* scale by alpha */          for (i = 0; i < nbx; i++) { /* scale by alpha */
535              REAL(BxP)[i] *= alpha;              REAL(BxP)[i] *= alpha;
# Line 549  Line 549 
549              BTp = Calloc(nrbB, int);              BTp = Calloc(nrbB, int);
550              triplet_to_col(ncbB, nrbB, nnz, tmp, Bi, Bx, BTp, BTi, BTx);              triplet_to_col(ncbB, nrbB, nnz, tmp, Bi, Bx, BTp, BTi, BTx);
551                                  /* sanity check */                                  /* sanity check */
552              if (BTp[nrbB] != nnz) error("cscb_trcbsm: transpose operation failed");              if (BTp[nrbB] != nnz) error(_("cscb_trcbsm: transpose operation failed"));
553              Free(tmp);              Free(tmp);
554                                  /* Solve one column at a time */                                  /* Solve one column at a time */
555              rhs = Calloc(ncbB, double);              rhs = Calloc(ncbB, double);
# Line 565  Line 565 
565                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);
566              }              }
567          }          }
568          error("cscb_trcbsm: method not yet written");          error(_("cscb_trcbsm: method not yet written"));
569      }      }
570      error("cscb_trcbsm: method not yet written");      error(_("cscb_trcbsm: method not yet written"));
571  }  }
572    
573  /**  /**
# Line 617  Line 617 
617          if (adims[1] != bdims[1] ||          if (adims[1] != bdims[1] ||
618              adims[0] != cdims[0] ||              adims[0] != cdims[0] ||
619              bdims[0] != cdims[1])              bdims[0] != cdims[1])
620              error("C[%d,%d,%d] := A[%d,%d,%d] %*% t(B[%d,%d,%d])",              error(_("C[%d,%d,%d] := A[%d,%d,%d] %*% t(B[%d,%d,%d])"),
621                    cdims[0], cdims[1], cdims[2],                    cdims[0], cdims[1], cdims[2],
622                    adims[0], adims[1], adims[2],                    adims[0], adims[1], adims[2],
623                    bdims[0], bdims[1], bdims[2]);                    bdims[0], bdims[1], bdims[2]);
624          if (nca != ncb)          if (nca != ncb)
625              error("C := A(ncblocks = %d) %*% t(B(ncblocks = %d)", nca, ncb);              error(_("C := A(ncblocks = %d) %*% t(B(ncblocks = %d)"), nca, ncb);
626          if (beta != 1.) {       /* scale C by beta */          if (beta != 1.) {       /* scale C by beta */
627              int ctot = cdims[0] * cdims[1] * cdims[2];              int ctot = cdims[0] * cdims[1] * cdims[2];
628              for (jj = 0; jj < ctot; jj++) Cx[jj] *= beta;              for (jj = 0; jj < ctot; jj++) Cx[jj] *= beta;
# Line 641  Line 641 
641          }          }
642          return;          return;
643      }      }
644      error("Code not yet written");      error(_("Code not yet written"));
645  }  }
646    
647  /**  /**

Legend:
Removed from v.565  
changed lines
  Added in v.582

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