kernelB.cl

#define NUM_NODES 25
#define NUM_ARCS  80
#define FLOAT_MAX 9999999
 
float setCost(uint arc, constant float *fftt, float *flow, constant float *cap)
{
  return fftt[arc] * (1.0F + 0.15F*pow((flow[arc]/cap[arc]), 4.0F));
}

float calcCost(uint arc,
               float deltaV,
               constant float *fftt,
               float *flow,
               constant float *cap)
{
  return fftt[arc] * (1.0F + 0.15F*pow(((flow[arc]+deltaV)/cap[arc]), 4.0F));
}

float calcCostPrime(uint arc, float deltaV, constant float *fftt, float *flow, constant float *cap)
{
  return 0.6F * fftt[arc] * pow(flow[arc] + deltaV, 3.0F) / pow(cap[arc], 4.0F);
}

int isInNodeArray(uint *array, uint n)
{
  uint i;
  for(i = 0; i < NUM_NODES; ++i) {
    if (array[i] == n) {
      return 1;
    }
  }
  return 0;
}

void updateCost(float *minCost,
                float *maxCost,
                float *cost,
                float *flow,
                uint sinkNode,
                constant uint *pointer,
                constant uint *tail,
                constant uint *head,
                uint *topoOrder,
                uint *nextMinArc,
                uint *nextMaxArc)
{
  uint i;
  uint arc;
  float costCand;

  for(i = 0; i < NUM_NODES; i++) {
    minCost[i] = FLOAT_MAX;
    maxCost[i] = 0;
  }

  minCost[sinkNode-1] = 0;
  nextMinArc[sinkNode-1] = NUM_ARCS;
  for(i = 0; i < NUM_NODES; i++) {
    arc = pointer[topoOrder[i]-1]-1;
    while(head[arc] == topoOrder[i]) {
      costCand = cost[arc] + minCost[head[arc]-1];
      if(costCand < minCost[tail[arc]-1]) {
        minCost[tail[arc]-1] = costCand;
        nextMinArc[tail[arc]-1] = arc;
      }
      ++arc;
    }
  }

  for(i = 0; i < NUM_NODES; i++) {
    maxCost[i] = minCost[i];
  }

  maxCost[sinkNode-1] = 0;
  nextMaxArc[sinkNode-1] = NUM_ARCS;
  for(i = 0; i < NUM_NODES; i++) {
    arc = pointer[topoOrder[i]-1]-1;
    while(head[arc] == topoOrder[i]) {
      if(cost[arc] >= FLOAT_MAX) {
        ++arc;
        continue;
      }
      costCand = cost[arc] + maxCost[head[arc]-1];
      if(costCand >= maxCost[tail[arc]-1] && flow[arc] > 0) {
        maxCost[tail[arc]-1] = costCand;
        nextMaxArc[tail[arc]-1] = arc;
      }
      ++arc;
    }
  }
}

void getPaths(uint origin,
              uint sinkNode,
              constant uint *head,
              uint *minPath,
              uint *maxPath,
              uint *nextMinArc,
              uint *nextMaxArc)
{
  uint i, j, k;
  uint length = 0;
  uint hasCommonNode = 1;

  minPath[0] = origin;
  maxPath[0] = origin;

  while(hasCommonNode) {
    if(minPath[0] == sinkNode - 1 || maxPath[0] == sinkNode - 1) {
      return;
    } else if(nextMaxArc[minPath[0]] == nextMinArc[maxPath[0]]) {
      minPath[0] = head[nextMaxArc[minPath[0]]]-1; 
      maxPath[0] = minPath[0]; 
      continue;
    }
    hasCommonNode = 0;
  }

  while(!hasCommonNode) {
    ++length;
    if(minPath[length - 1] != sinkNode-1 && minPath[length - 1] != NUM_NODES) {
      minPath[length] = head[nextMinArc[minPath[length-1]]]-1;
    } else {
      minPath[length] = NUM_NODES;
    }
    if(maxPath[length - 1] != sinkNode-1 && maxPath[length - 1] != NUM_NODES) {
      maxPath[length] = head[nextMaxArc[maxPath[length-1]]]-1;
    } else {
      maxPath[length] = NUM_NODES;
    }
    for(i = 1; i < length + 1; ++i) {
      for(j = 1; j < length + 1; ++j) {
        if(maxPath[i] == minPath[j]) {
          for(k = j + 1; k < NUM_NODES; ++k) {
            minPath[k] = NUM_NODES;
          }
          for(k = i + 1; k < NUM_NODES; ++k) {
            maxPath[k] = NUM_NODES;
          }
          hasCommonNode = 1;
          return;
        }
      }
    }
  }
}

void shiftFlow(uint sinkNode,
               uint *maxPath,
               uint *minPath,
               constant uint *pointer,
               constant uint *tail,
               constant uint *head,
               uint *topoOrder,
               uint *nextMinArc,
               uint *nextMaxArc,
               float *flow,
               constant float *fftt,
               constant float *cap,
               float *cost,
               float *minCost,
               float *maxCost)
{
  uint i = 0;
  float nextFlow;
  float max = FLOAT_MAX;
  float root = 0;
  float f = 1.0F;
  float fPrime;
  uint maxIter = 10;
  uint iter = 0;

  for(i = 0; i < NUM_NODES; ++i) {
    if(nextMaxArc[maxPath[i]] >= NUM_ARCS
      || maxPath[i] >= NUM_NODES) {
        break;
    }
    nextFlow = flow[nextMaxArc[maxPath[i]]];
    if(nextFlow < max) {
      max = nextFlow;
    }
  }
  root = max;

  while(fabs(f) > 0.0001F) {
    if(iter >= maxIter) {
      break;
    }
    ++iter;
    i = 0;
    f = 0; 
    fPrime = 0; 
    while(minPath[i+1] < NUM_NODES || maxPath[i+1] < NUM_NODES) {
      if(minPath[i+1] < NUM_NODES) {
        f += calcCost(nextMinArc[minPath[i]], root, fftt, flow, cap);
        fPrime += calcCostPrime(nextMinArc[minPath[i]], root, fftt, flow, cap);
      }
      if(maxPath[i+1] < NUM_NODES) {
        f -= calcCost(nextMaxArc[maxPath[i]], -1*root, fftt, flow, cap);
        fPrime += calcCostPrime(nextMaxArc[maxPath[i]], -1*root, fftt, flow, cap);
      }
      ++i;
    }
    root -= f/fPrime;
  }
  if(root > max) {
    root = max;
  } else if(root < 0.0001f) {
    root = 0.1 * max;
    //return;
  } 
  i = 0;
  while(minPath[i+1] < NUM_NODES || maxPath[i+1] < NUM_NODES) {
    if(minPath[i+1] < NUM_NODES) {
      flow[nextMinArc[minPath[i]]] += root;
      cost[nextMinArc[minPath[i]]] = setCost(nextMinArc[minPath[i]], fftt, flow, cap);
    }
    if(maxPath[i+1] < NUM_NODES) {
      flow[nextMaxArc[maxPath[i]]] -= root;
      cost[nextMaxArc[maxPath[i]]] = setCost(nextMaxArc[maxPath[i]], fftt, flow, cap);
    }
    ++i;
  }
  updateCost(minCost, maxCost, cost, flow, sinkNode, pointer, tail, head, topoOrder, nextMinArc, nextMaxArc);
}

void equilibrateBush(uint sinkNode, 
                     uint *maxPath,
                     uint *minPath,
                     uint *nextMinArc,
                     uint *nextMaxArc,
                     constant uint *pointer,
                     float *demand,
                     constant uint *tail,
                     constant uint *head,
                     uint *topoOrder,
                     float *flow,
                     constant float *fftt,
                     constant float *cap,
                     float *cost,
                     float *minCost,
                     float *maxCost)
{
  int i;
  int j;
  int eqBush = 0;
  int updated;
  int iter = 0;
  float tol = 0.001f;

  while(!eqBush) {
    updated = 0;
    if(iter > 0 && iter % 100 == 0) {
      tol = tol * 10.f;
    }
    for(i = 0; i < NUM_NODES; ++i) {
      j = topoOrder[i] - 1;
      if(fabs(maxCost[j] - minCost[j]) > tol && demand[j] > 0) {
        getPaths(j, sinkNode, head, minPath, maxPath, nextMinArc, nextMaxArc);
        shiftFlow(sinkNode, maxPath, minPath, pointer, tail, head, topoOrder, nextMinArc, nextMaxArc, flow, fftt, cap, cost, minCost, maxCost);
        updated = 1;
      }
    }
    if(updated == 0) {
      eqBush = 1;
    }
    iter++;
  }
}

void initNetwork(const uint sinkNode,
                 uint *topoOrder,
                 constant uint *tail,
                 constant uint *head,
                 constant uint *pointer,
                 float *demand,
                 uint *nextMinArc,
                 uint *nextMaxArc,
                 float *minCost,
                 float *maxCost,
                 float *cost,
                 constant float *fftt,
                 constant float *cap,
                 float *flow)
{
  uint i;
  uint place;
  uint queue;
  uint arc;
  float costCand;
  uint nextArc;

  // set topological order of nodes
  for(i = 0; i < NUM_NODES; ++i) {
    topoOrder[i] = 0;
  }
  topoOrder[0] = sinkNode;
  queue = 1;
  for(place = 0; place < NUM_NODES; ++place) {
    for(i=0; i < NUM_ARCS; ++i) {
      if (head[i] == topoOrder[place] && !isInNodeArray(topoOrder, tail[i])) {
        topoOrder[queue] = tail[i];
        ++queue;
      }
    }
  }

  // init cost
  for(i = 0; i < NUM_ARCS; ++i) {
    cost[i] = FLOAT_MAX;
  }

  // init minimum paths
  for(i = 0; i < NUM_NODES; i++) {
    minCost[i] = FLOAT_MAX;
  }
  minCost[sinkNode-1] = 0;
  nextMinArc[sinkNode-1] = NUM_ARCS;
  for(i = 0; i < NUM_NODES; i++) {
    arc = pointer[topoOrder[i]-1]-1;
    while(head[arc] == topoOrder[i]) {
      costCand = fftt[arc] + minCost[head[arc]-1];
      if(costCand < minCost[tail[arc]-1]) {
        minCost[tail[arc]-1] = costCand;
        nextMinArc[tail[arc]-1] = arc;
      }
      ++arc;
    }
  }

  // init maximum paths
  for(i = 0; i < NUM_NODES; i++) {
    maxCost[i] = minCost[i];
  }

  // init flows on bush
  nextArc = NUM_ARCS;
  for(i = 0; i < NUM_NODES; i++) {
    nextArc = nextMinArc[i];
    while(nextArc < NUM_ARCS) {
      flow[nextArc] += demand[i];
      nextArc = nextMinArc[head[nextArc]-1];
    }
  }

  // init costs on bush
  for(arc = 0; arc < NUM_ARCS; ++arc) {
    if(minCost[head[arc]-1] < minCost[tail[arc]-1]) {
      cost[arc] = setCost(arc, fftt, flow, cap);
    }
  }

  updateCost(minCost, maxCost, cost, flow, sinkNode, pointer, tail, head, topoOrder, nextMinArc, nextMaxArc);
}

int updateBush(uint sinkNode,
               uint topoOrder[],
               constant uint tail[],
               constant uint head[],
               constant uint pointer[],
               uint nextMinArc[],
               uint nextMaxArc[],
               float minCost[],
               float maxCost[],
               float cost[],
               constant float *fftt,
               constant float *cap,
               float flow[])

{
  uint i;
  uint j;
  uint isUpdated = 0;

  for(i = 0; i < NUM_ARCS; ++i) {
    if(flow[i] == 0 && cost[i] < FLOAT_MAX) {
      j = pointer[tail[i]-1]-1;
      while(head[j] == tail[i]) {
        if(tail[j] == head[i] && fftt[j] + minCost[head[j]-1] < minCost[tail[j]-1]) {
          cost[i] = FLOAT_MAX;
          cost[j] = setCost(j, fftt, flow, cap); 
          isUpdated = 1;
        }
        ++j;
      }
    }
  }
  updateCost(minCost, maxCost, cost, flow, sinkNode, pointer, tail, head, topoOrder, nextMinArc, nextMaxArc);
  return isUpdated;
}

void equilibrateNetwork(uint sinkNode,
                        uint *maxPath,
                        uint *minPath,
                        uint *nextMinArc,
                        uint *nextMaxArc,
                        constant uint *pointer,
                        float *demand,
                        constant uint *tail,
                        constant uint *head,
                        uint *topoOrder,
                        float *flow,
                        constant float *fftt,
                        constant float *cap,
                        float *cost,
                        float *minCost,
                        float *maxCost)
{
  uint bushUpdated = 1;
  while(bushUpdated) {
    equilibrateBush(sinkNode, maxPath, minPath, nextMinArc, nextMaxArc, pointer, demand, tail, head, topoOrder, flow, fftt, cap, cost, minCost, maxCost);
    bushUpdated = updateBush(sinkNode, topoOrder, tail, head, pointer, nextMinArc, nextMaxArc, minCost, maxCost, cost, fftt, cap, flow); 
  }
}

__kernel void kernelB(__global unsigned int * output,
                      __global unsigned int * input,
                      constant unsigned int * pointer,
                      constant unsigned int * tail,
                      constant unsigned int * head,
                      constant float * fftt,
                      constant float * cap)

{

  uint sinkNode;
  uint tid = get_global_id(0);
  sinkNode = (tid % 25) + 1;

  // float demand[] = {50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 50, 0, 50};
  // float demand[] = {300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 300, 0, 300, 300, 300, 300, 300, 300, 300, 300, 0, 300};
  float demand[] = {200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200, 200};


  // arc constants

  // node variables
  uint topoOrder[NUM_NODES];
  float minCost[NUM_NODES];
  float maxCost[NUM_NODES];
  uint nextMaxArc[NUM_NODES];
  uint nextMinArc[NUM_NODES];
  uint minPath[NUM_NODES];
  uint maxPath[NUM_NODES];

  // arc variables
  float flow[NUM_ARCS];
  float cost[NUM_ARCS];

  initNetwork(sinkNode, topoOrder, tail, head, pointer, demand, nextMinArc, nextMaxArc, minCost, maxCost, cost, fftt, cap, flow);
  equilibrateNetwork(sinkNode, maxPath, minPath, nextMinArc, nextMaxArc, pointer, demand, tail, head, topoOrder, flow, fftt, cap, cost, minCost, maxCost);

  output[tid] = (int) maxCost[0];
}