spsm.cc

sparselib库

    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_UpperDiagSolve_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] + 1; j < pntr[i+1]; j++)
      z += c[indx[j]] * val[j];
    c[i+pntr[0]] = (b[i] - z) / val[pntr[i]];
      }
    else
      for (l = 0; l < n; l++)
    for (i = m - 1; i >= 0; i--) {  
      z = 0;
      for (j = pntr[i] + 1; 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) / val[pntr[i]];
    }
  } else {
    if (n == 1)
      for (i = m - 1; i >= 0; i--) {    
    z = 0;
    for (j = pntr[i] + 1; j < pntr[i+1]; j++)
      z += c[indx[j]] * val[j];
    c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
      }
    else
      for (l = 0; l < n; l++)
    for (i = m - 1; i >= 0; i--) {  
      z = 0;
      for (j = pntr[i] + 1; 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) / val[pntr[i]];
    }
  }
}


static void
CompRow_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;

  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] + 1; j < pntr[i+1]; j++)
      z += c[indx[j]] * val[j];
    c[i+pntr[0]] = (b[i] - z) / val[pntr[i]];
      }
    else
      for (l = 0; l < n; l++)
    for (i = m - 1; i >= 0; i--) {  
      z = 0;
      for (j = pntr[i] + 1; 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) / val[pntr[i]];
    }
  } else {
    if (n == 1)
      for (i = m - 1; i >= 0; i--) {    
    z = 0;
    for (j = pntr[i] + 1; j < pntr[i+1]; j++)
      z += c[indx[j]] * val[j];
    c[i+pntr[0]] = alpha * (b[i] - z) / val[pntr[i]];
      }
    else
      for (l = 0; l < n; l++)
    for (i = m - 1; i >= 0; i--) {  
      z = 0;
      for (j = pntr[i] + 1; 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) / val[pntr[i]];
    }
  }
}







/*
 * C <- alpha D A^{-1} B + beta C
 * C <- alpha A^{-1} D B + beta C
 */
void F77NAME(scscsm)
  (const int &transa, const int &m, const int &n,
   const int &unitd, const float *dv, 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);
  }

  if (beta != 0.0) {
    std::cerr << "beta != 0 not implemented" << "\n";
    exit(1);
  }

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

  if (transa == 0) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    }
  } else if (transa == 1 || transa == 2) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else {
      std::cerr << "Bad descra[1]" << "\n";
      exit(1);
    }
  } else {
    std::cerr << "Bad transa" << "\n";
    exit(1);
  }
}


void F77NAME(scsrsm)
  (const int &transa, const int &m, const int &n,
   const int &unitd, const float *dv, 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);
  }

  if (beta != 0.0) {
    std::cerr << "beta != 0 not implemented" << "\n";
    exit(1);
  }

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

  if (transa == 0) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompRow_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompRow_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    }
  } else if (transa == 1 || transa == 2) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompCol_UpperDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_UpperUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompCol_LowerDiagSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_LowerUnitSolve_float(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else {
      std::cerr << "Bad descra[1]" << "\n";
      exit(1);
    }
  } else {
    std::cerr << "Bad transa" << "\n";
    exit(1);
  }
}


void F77NAME(dcscsm)
  (const int &transa, const int &m, const int &n,
   const int &unitd, const double *dv, 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);
  }

  if (beta != 0.0) {
    std::cerr << "beta != 0 not implemented" << "\n";
    exit(1);
  }

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

  if (transa == 0) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    }
  } else if (transa == 1 || transa == 2) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else {
      std::cerr << "Bad descra[1]" << "\n";
      exit(1);
    }
  } else {
    std::cerr << "Bad transa" << "\n";
    exit(1);
  }
}




void F77NAME(dcsrsm)
  (const int &transa, const int &m, const int &n,
   const int &unitd, const double *dv, 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);
  }

  if (beta != 0.0) {
    std::cerr << "beta != 0 not implemented" << "\n";
    exit(1);
  }

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

  if (transa == 0) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompRow_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompRow_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompRow_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    }
  } else if (transa == 1 || transa == 2) {
    if (descra[1] == 1) {
      if (descra[2] == 0)
    CompCol_UpperDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_UpperUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else if (descra[1] == 2) {
      if (descra[2] == 0)
    CompCol_LowerDiagSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else if (descra[2] == 1)
    CompCol_LowerUnitSolve_double(m, n, unitd, dv, alpha, val, indx, pntr,
                   b, ldb, c, ldc);
      else {
    std::cerr << "Bad descra[2]" << "\n";
    exit(1);
      }
    } else {
      std::cerr << "Bad descra[1]" << "\n";
      exit(1);
    }
  } else {
    std::cerr << "Bad transa" << "\n";
    exit(1);
  }
}