# SCM Repository

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

# Diff of /pkg/src/dgBCMatrix.c

revision 556, Wed Feb 16 20:49:44 2005 UTC revision 559, Fri Feb 18 14:29:03 2005 UTC
# Line 35  Line 35
35   * where A is a compressed, sparse, blocked matrix and   * where A is a compressed, sparse, blocked matrix and
36   * B and C are dense matrices.   * B and C are dense matrices.
37   *   *
38   * @param side 'L' or 'l' for left, otherwise right   * @param side LFT or RGT
39   * @param transa 'T' or 't' for transpose.   * @param transa TRN or NTR
40   * @param m number of rows in C   * @param m number of rows in C
41   * @param n number of columns in C   * @param n number of columns in C
42   * @param k number of rows in B if side = 'L', otherwise   * @param k number of rows in B if side == LFT, otherwise
43   *        number of columns in B.   *        number of columns in B.
44   * @param alpha   * @param alpha
45   * @param A pointer to a dgBCMatrix object   * @param A pointer to a dgBCMatrix object
# 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 == 'L'");          error("Call to cscb_mm must have side == LFT");
115      }      }
116  }  }
117
# Line 373  Line 373
373   * Solve a triangular system of the form op(A)*X = alpha*B where A   * Solve a triangular system of the form op(A)*X = alpha*B where A
374   * is a dgBCMatrix triangular matrix and B is a dense matrix.   * is a dgBCMatrix triangular matrix and B is a dense matrix.
375   *   *
376   * @param uplo 'U' or 'L' for upper or lower   * @param uplo UPP or LOW
377   * @param trans 'T' or 'N' for transpose or no transpose   * @param trans TRN or NTR
378   * @param diag 'U' or 'N' for unit diagonal or non-unit   * @param diag UNT or NUN
379   * @param A pointer to a triangular dgBCMatrix object   * @param A pointer to a triangular dgBCMatrix object
* @param B pointer to the contents of the matrix B
380   * @param m number of rows in B   * @param m number of rows in B
381   * @param n number of columns in B   * @param n number of columns in B
382     * @param B pointer to the contents of the matrix B
383   * @param ldb leading dimension of B as declared in the calling function   * @param ldb leading dimension of B as declared in the calling function
384   */   */
385  void  void
386  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa,  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa,
387            enum CBLAS_DIAG diag, double alpha, SEXP A,            enum CBLAS_DIAG diag, double alpha, SEXP A,
388            double B[], int m, int n, int ldb)            int m, int n, double B[], int ldb)
389  {  {
390      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
391          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
# Line 422  Line 422
422              }              }
423              return;              return;
424          } else {          } else {
425              int p, p2, sza = xdims[0] * xdims[0], szb = xdims[0] * n;              int p, p2, sza = xdims[0] * xdims[0];
426              double *tmp = Calloc(szb, double);
427              if (uplo == UPP) error("Code for upper triangle not yet written");              if (uplo == UPP) error("Code for upper triangle not yet written");
428              if (transa == TRN) {              if (transa == TRN) {
429                  for (j = nb - 1; j >= 0; j--) {                  for (j = nb - 1; j >= 0; j--) {
430                      p2 = Ap[j+1];                      p2 = Ap[j+1];

F77_CALL(dlacpy)("A", xdims, &n, B + j * xdims[0], &ldb,
tmp, xdims);
431                      for (p = Ap[j]; p < p2; p++)                      for (p = Ap[j]; p < p2; p++)
432                          F77_CALL(dgemm)("T", "N", xdims, &n, xdims,                          F77_CALL(dgemm)("T", "N", xdims, &n, xdims,
433                                          &minus1, Ax + p * sza, xdims,                                          &minus1, Ax + p * sza, xdims,
434                                          B + Ai[p] * xdims[0], &ldb,                                          B + Ai[p] * xdims[0], &ldb,
435                                          &one, tmp, xdims);                                          &one, B + j * xdims[0], &ldb);
F77_CALL(dlacpy)("A", xdims, &n, tmp, xdims,
B + j * xdims[0], &ldb);
436                  }                  }
437              } else {              } else {
438                  for (j = 0; j < nb; j++) {                  for (j = 0; j < nb; j++) {
439                      p2 = Ap[j+1];                      p2 = Ap[j+1];

F77_CALL(dlacpy)("A", xdims, &n, B + j * xdims[0], &ldb,
tmp, xdims);
440                      for (p = Ap[j]; p < p2; p++)                      for (p = Ap[j]; p < p2; p++)
441                          F77_CALL(dgemm)("N", "N", xdims, &n, xdims,                          F77_CALL(dgemm)("N", "N", xdims, &n, xdims,
442                                          &minus1, Ax + p * sza, xdims,                                          &minus1, Ax + p * sza, xdims,
443                                          B + Ai[p] * xdims[0], &ldb,                                          B + j * xdims[0], &ldb,
444                                          &one, tmp, xdims);                                          &one, B + Ai[p] * xdims[0], &ldb);
F77_CALL(dlacpy)("A", xdims, &n, tmp, xdims,
B + j * xdims[0], &ldb);
445                  }                  }
446              }              }
447          }          }

Legend:
 Removed from v.556 changed lines Added in v.559