spmm.cc

http://gams.cam.nist.gov/acmd/Staff/RPozo/sparselib++.html

      for (i = 0; i < n; i++)
    for (j = 0; j < m; j++)
      _SpMatVal(c, ldc, j, i) *= beta;
  }
}

/*
 * dcoom -- coordinate format matrix-matrix multiply
 *
 * C <- alpha A B + beta C
 *
 * Arguments:
 *
 * int &transa  Indicates how to operate with the sparse matrix
 *      0 : operate with matrix
 *      1 : operate with transpose matrix
 *      2 : operate with conjugate transpose matrix
 *
 * int &m   Number of rows in matrix c
 *
 * int &n   Number of columns in matrix c
 *
 * int &k   Number of rows in matrix b
 *
 * double &alpha Scalar parameter
 *
 * double &beta Scalar parameter
 *
 * int descra[] Descriptor argument.  Nine element integer array
 *      descra[0] matrix structure
 *          0 : general
 *          1 : symmetric
 *          2 : Hermition
 *          3 : Triangular
 *          4 : Anti-Symmetric
 *          5 : Diagonal
 *      descra[1] upper/lower triangular indicator
 *          1 : lower
 *          2 : upper
 *      descra[2] main diagonal type
 *          0 : non-unit
 *          1 : unit
 *      descra[4] repeated indices?
 *          0 : unknown
 *          1 : no repeated indices
 *
 *
 * double *val  scalar array of length nnz containing matrix entries
 *
 * int *indx    integer array of length nnz containing row indices
 *
 * int *jndx    integer array of length nnz containing column indices
 *
 * double *b    rectangular array with first dimension ldb
 *
 * double *c    rectangular array with first dimension ldc
 *
 * double *work scratch array of length lwork.  lwork should be at least
 *      max(m,n)
 *
 */
void F77NAME(scoomm)
  (const int &transa, const int &m, const int &n, const int &k,
   const float &alpha,
   const int descra[], const float *val,
   const int *indx, const int *jndx, const int &nnz,
   const float *b, const int &ldb,
   const float &beta, float *c, const int &ldc,
   float *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_float(m, n, c, ldc, beta);

  if (alpha == 0.0)
    return;

  // Use this hack if transpose is desired
  if (transa == 1 || transa == 2) {
    const int *itmp = indx;
    indx = jndx;
    jndx = itmp;
  }
  CoordMatVec_float(m, n, k, alpha, val, indx, jndx, nnz, b, ldb, c, ldc);
}


void F77NAME(dcoomm)
  (const int &transa, const int &m, const int &n, const int &k,
   const double &alpha,
   const int descra[], const double *val,
   const int *indx, const int *jndx, const int &nnz,
   const double *b, const int &ldb,
   const double &beta, double *c, const int &ldc,
   double *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_double(m, n, c, ldc, beta);

  if (alpha == 0.0)
    return;

  // Use this hack if transpose is desired
  if (transa == 1 || transa == 2) {
    const int *itmp = indx;
    indx = jndx;
    jndx = itmp;
  }
  CoordMatVec_double(m, n, k, alpha, val, indx, jndx, nnz, b, ldb, c, ldc);
}




/*
 * dcscm -- comp sparse column matrix-matrix multiply
 *
 * Arguments:
 *
 * int &transa  Indicates how to operate with the sparse matrix
 *      0 : operate with matrix
 *      1 : operate with transpose matrix
 *      2 : operate with conjugate transpose matrix
 *
 * int &m   Number of rows in matrix c
 *
 * int &n   Number of columns in matrix c
 *
 * int &k   Number of rows in matrix b
 *
 * double &alpha Scalar parameter
 *
 * double &beta Scalar parameter
 *
 * int descra[] Descriptor argument.  Nine element integer array
 *      descra[0] matrix structure
 *          0 : general
 *          1 : symmetric
 *          2 : Hermition
 *          3 : Triangular
 *          4 : Anti-Symmetric
 *          5 : Diagonal
 *      descra[1] upper/lower triangular indicator
 *          1 : lower
 *          2 : upper
 *      descra[2] main diagonal type
 *          0 : non-unit
 *          1 : unit
 *
 * double *val  scalar array of length nnz containing matrix entries
 *
 * int *indx    integer array of length nnz containing row indices
 *
 * int *pntr    integer array of length k+1 such that pntr(j)-pntr(1)
 *      points to location in val of the first element in column j
 *
 * double *b    rectangular array with first dimension ldb
 *
 * double *c    rectangular array with first dimension ldc
 *
 * double *work scratch array of length lwork.  lwork should be at least
 *      max(m,n)
 *
 */
void F77NAME(scscmm)
  (const int &transa, const int &m, const int &n, const int &k,
   const float &alpha,
   const int descra[], const float *val,
   const int *indx, const int *pntr, const float *b, int &ldb,
   const float &beta, float *c, const int &ldc,
   float *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_float(m, n, c, ldc, beta);

  if (transa == 1 || transa == 2)
    CompRowMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
  else
    CompColMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
}


void F77NAME(dcscmm)
  (const int &transa, const int &m, const int &n, const int &k,
   const double &alpha,
   const int descra[], const double *val,
   const int *indx, const int *pntr, const double *b, int &ldb,
   const double &beta, double *c, const int &ldc,
   double *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_double(m, n, c, ldc, beta);

  if (transa == 1 || transa == 2)
    CompRowMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
  else
    CompColMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
}




/*
 * dcsrm -- comp sparse row matrix-matrix multiply
 *
 * Arguments:
 *
 * int &transa  Indicates how to operate with the sparse matrix
 *      0 : operate with matrix
 *      1 : operate with transpose matrix
 *      2 : operate with conjugate transpose matrix
 *
 * int &m   Number of rows in matrix c
 *
 * int &n   Number of columns in matrix c
 *
 * int &k   Number of rows in matrix b
 *
 * double &alpha Scalar parameter
 *
 * double &beta Scalar parameter
 *
 * int descra[] Descriptor argument.  Nine element integer array
 *      descra[0] matrix structure
 *          0 : general
 *          1 : symmetric
 *          2 : Hermition
 *          3 : Triangular
 *          4 : Anti-Symmetric
 *          5 : Diagonal
 *      descra[1] upper/lower triangular indicator
 *          1 : lower
 *          2 : upper
 *      descra[2] main diagonal type
 *          0 : non-unit
 *          1 : unit
 *
 * double *val  scalar array of length nnz containing matrix entries
 *
 * int *indx    integer array of length nnz containing column indices
 *
 * int *pntr    integer array of length k+1 such that pntr(j)-pntr(1)
 *      points to location in val of the first element in row j
 *
 * double *b    rectangular array with first dimension ldb
 *
 * double *c    rectangular array with first dimension ldc
 *
 * double *work scratch array of length lwork.  lwork should be at least
 *      max(m,n)
 *
 */
void F77NAME(scsrmm)
  (const int &transa, const int &m, const int &n, const int &k,
   const float &alpha,
   const int descra[], const float *val,
   const int *indx, const int *pntr, const float *b, int &ldb,
   const float &beta, float *c, const int &ldc,
   float *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_float(m, n, c, ldc, beta);

  if (transa == 1 || transa == 2)
    CompColMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
  else
    CompRowMatVec_float(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
}


void F77NAME(dcsrmm)
  (const int &transa, const int &m, const int &n, const int &k,
   const double &alpha,
   const int descra[], const double *val,
   const int *indx, const int *pntr, const double *b, int &ldb,
   const double &beta, double *c, const int &ldc,
   double *work, const int &lwork)
{
  if (descra[0] != 0) {
    std::cerr << "Must have general matrix" << "\n";
    exit(1);
  }

  // To make the compiler happy
  if (work && lwork)
    ;

  ScaleRectangularArray_double(m, n, c, ldc, beta);

  if (transa == 1 || transa == 2)
    CompColMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
  else
    CompRowMatVec_double(m, n, k, alpha, val, indx, pntr, b, ldb, c, ldc);
}