spsm.cc

sparselib库

    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 void
CompCol_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 void
CompRow_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 void
CompRow_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 void
CompRow_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 void
CompRow_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 void
CompRow_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 void
CompRow_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);
      }