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);}