std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_double(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = c[i+pntr[0]] / val[pntr[i+1]-1];    c[i+pntr[0]] = z;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = alpha * c[i+pntr[0]] / val[pntr[i+1]-1];    c[i+pntr[0]] = z;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompCol_UpperDiagSolve_float(int m, int n, int unitd, const float *dv,                 float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  CopyRectangularArray_float(m, n, b, ldb, &c[pntr[0]-1], ldc);  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = c[i+pntr[0]] / val[pntr[i+1]-1];    c[i+pntr[0]] = z;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = alpha * c[i+pntr[0]] / val[pntr[i+1]-1];    c[i+pntr[0]] = z;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      c[indx[j]] -= z * val[j];      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = alpha * _SpMatVal(c, ldc, i+pntr[0], l) / val[pntr[i+1]-1];      _SpMatVal(c, ldc, i+pntr[0], l) = z;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        _SpMatVal(c, ldc, indx[j], l) -= z * val[j];    }  }}static voidCompRow_LowerUnitSolve_double(int m, int n, int unitd, const double *dv,                 double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = b[i] - z;      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+1] = alpha * (b[i] - z);      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);    }  }}static voidCompRow_LowerUnitSolve_float(int m, int n, int unitd, const float *dv,                 float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = b[i] - z;      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+1] = alpha * (b[i] - z);      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);    }  }}static voidCompRow_LowerDiagSolve_double(int m, int n, int unitd, const double *dv,                 double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = (b[i] - z) / val[pntr[i+1]-1];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) =         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i+1]-1];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) =         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];    }  }}static voidCompRow_LowerDiagSolve_float(int m, int n, int unitd, const float *dv,                 float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = (b[i] - z) / val[pntr[i+1]-1];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) =         (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];    }  } else {    if (n == 1)      for (i = 0; i < m; i++) {    z = 0;    for (j = pntr[i]; j < pntr[i+1] - 1; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i+1]-1];      }    else      for (l = 0; l < n; l++)    for (i = 0; i < m; i++) {      z = 0;      for (j = pntr[i]; j < pntr[i+1] - 1; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) =         alpha * (_SpMatVal(b, ldb, i, l) - z) / val[pntr[i+1]-1];    }  }}static voidCompRow_UpperUnitSolve_double(int m, int n, int unitd, const double *dv, double alpha,               const double *val, const int *indx, const int *pntr,               const double *b, int ldb, double *c, int ldc){  int i, j, l;  double z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = b[i] - z;      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = alpha * (b[i] - z);      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = alpha * (_SpMatVal(b, ldb, i, l) - z);    }  }}static voidCompRow_UpperUnitSolve_float(int m, int n, int unitd, const float *dv, float alpha,               const float *val, const int *indx, const int *pntr,               const float *b, int ldb, float *c, int ldc){  int i, j, l;  float z;  // To make the compiler happy  if (dv)    ;  if (unitd != 1) {    std::cerr << "unitd != 1 not implemented" << "\n";    exit(1);  }  if (alpha == 0.0)    return;  c -= 1;  val -= pntr[0];  indx -= pntr[0];  if (alpha == 1.0) {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = b[i] - z;      }    else      for (l = 0; l < n; l++)    for (i = m - 1; i >= 0; i--) {      z = 0;      for (j = pntr[i]; j < pntr[i+1]; j++)        z += _SpMatVal(c, ldc, indx[j], l) * val[j];      _SpMatVal(c, ldc, i+pntr[0], l) = (_SpMatVal(b, ldb, i, l) - z);    }  } else {    if (n == 1)      for (i = m - 1; i >= 0; i--) {    z = 0;    for (j = pntr[i]; j < pntr[i+1]; j++)      z += c[indx[j]] * val[j];    c[i+pntr[0]] = alpha * (b[i] - z);      }