diff --git a/.gitignore b/.gitignore index 87236f9..26adc12 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +msvc build Debug Release diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index ce17782..e51990e 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -84,44 +84,17 @@ include_directories( ${METIS_PATH}/Lib ) -set(CORE_CUDA_SOURCES - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/aggregator.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/mis.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/randMIS.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/amg.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/amg_config.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/amg_level.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/amg_signal.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/cgcycle.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/allocator.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/ComputePermutationMethods.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/cutil.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/smoothedMG_amg_level.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/smoother.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/gauss_seidel.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/FEM2D.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/FEM3D.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/tetmesh.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/AggMIS_Types.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/misHelpers.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/randomizedMIS_GPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/setup_solver.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/cuda_resources.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_Aggregation_CPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_Aggregation_GPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_GraphHelpers.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_IOHelpers.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_MergeSplitConditioner_CPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_MergeSplitConditioner.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_MIS_CPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/AggMIS_MIS_GPU.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/TriMesh_connectivity.cu" - "${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/TriMesh_io.cu" -) - -set(CORE_CUDA_HEADERS - "${CMAKE_CURRENT_SOURCE_DIR}/cuda/perform_element_loop_3D.cuh" -) +FILE(GLOB CORE_CUDA_SOURCES ${CMAKE_CURRENT_SOURCE_DIR}/cuda/*.cu + ${CMAKE_CURRENT_SOURCE_DIR}/aggmis/cuda/*.cu) + +FILE(GLOB CORE_CUDA_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/include/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/smoothers/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/smoothedMG/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/smoothedMG/aggregators/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/FEM/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/include/cycles/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/aggmis/include/*.h + ${CMAKE_CURRENT_SOURCE_DIR}/cuda/perform_element_loop_3D.cuh) CUDA_ADD_LIBRARY(FEM_CORE STATIC ${CORE_CUDA_SOURCES} ${CORE_CUDA_HEADERS}) ADD_DEPENDENCIES(FEM_CORE metis) diff --git a/src/core/aggmis/cuda/AggMIS_Testing.cu b/src/core/aggmis/cuda/AggMIS_Testing.cu deleted file mode 100644 index 71965c2..0000000 --- a/src/core/aggmis/cuda/AggMIS_Testing.cu +++ /dev/null @@ -1,665 +0,0 @@ -#include "AggMIS_Testing.h" -namespace AggMIS { - namespace Testing { - MetricsContext::MetricsContext(Graph_CSR& graph, IdxVector_h& aggregation) { - this->graph = &graph; - this->aggregation.resize(aggregation.size()); - thrust::copy(aggregation.begin(), aggregation.end(), this->aggregation.begin()); - Initialize(); - } - MetricsContext::MetricsContext(int size, int* adjIndices, int* adjacency, int* aggregation) { - graph = new Graph_CSR(size, adjIndices, adjacency); - Initialize(); - } - - // Public Methods - double MetricsContext::GetConvexityRatio(int aggId) { - if (aggId < aggregates.size() && aggId >= 0) - { - distanceLookups = 0; - makeConvexCalls = 0; - EnsureConvex(aggId); - printf("Found convexity ratio for aggregate %d of size %d with %d MakeConvex calls and %d FindDistances calls.\n", aggId, aggregates[aggId].size(), makeConvexCalls, distanceLookups); - - return (double)aggregates[aggId].size() / convexAggregates[aggId].size(); - } - return -1.0; - } - double MetricsContext::GetEccentricityRatio(int aggId) { - if (aggId < aggregates.size() && aggId >= 0) - { - EnsureConvex(aggId); - return GetEccentricityRatio(convexAggregates[aggId]); - } - return -1.0; - } - int MetricsContext::GetAggregateCount() { - return aggregates.size(); - } - - // Internal Methods - double MetricsContext::GetEccentricityRatio(vector& aggregate) { - // Find the centroidal nodes - vector *centroid = FindCentroid(aggregate); - - // Find the distances from the centroidal nodes - vector*> distances(centroid->size()); - for (int i = 0; i < distances.size(); i++) - distances[i] = FindDistances(centroid->at(i), aggregate); - - // Find farthest internal node by finding highest average distance to centroidal nodes - double farthestInternal = 0.0; - for (int i = 0; i < aggregate.size(); i++) - { - double distance = 0.0; - for (int j = 0; j < distances.size(); j++) - distance += distances[j]->at(aggregate[i]); - distance /= distances.size(); - farthestInternal = std::max(distance, farthestInternal); - } - - // Find external nodes seen by centroidal nodes - set externalCandidates; - for (int i = 0; i < distances.size(); i++) - { - map::iterator it; - for (it = distances[i]->begin(); it != distances[i]->end(); it++) - { - int node = (*it).first; - if(!std::binary_search(aggregate.begin(), aggregate.end(), (*it).first)) - externalCandidates.insert(node); - } - } - - // Check each external to find the one with minimum average distance - double nearestExternal = 100000; - set::iterator sIt; - for (sIt = externalCandidates.begin(); sIt != externalCandidates.end(); sIt++) - { - double currentScore = 0.0; - bool unseen = false; - for (int i = 0; i < distances.size(); i++) - { - if (distances[i]->count(*sIt) > 0) - currentScore += distances[i]->at(*sIt); - else - unseen = true; - } - if (!unseen) - { - currentScore /= distances.size(); - if (currentScore < nearestExternal) - nearestExternal = currentScore; - } - } - - // Memory clean up: - for (int i = 0; i < distances.size(); i++) - delete distances[i]; - - // Return value for metric - return (nearestExternal - 1) / farthestInternal; - } - double MetricsContext::GetMinimumEnclosingBallRatio(vector& aggregate) { - /* Method: find distances out from first and second aggregate nodes then - find the intersection of the distance envelopes, then find the next nodes - envelope and intersect it with the first intersection. The center of the - minimal enclosing ball must lie within the intersection of all nodes distance - envelopes*/ - - // Ensure input it sorted - std::sort(aggregate.begin(), aggregate.end()); - - // If there is only one node in aggregate return - if (aggregate.size() < 2) - return 1.0; - - // Otherwise find the intersection of the first two nodes distance envelope - map *a = FindDistances(aggregate[0], aggregate); - map *b = FindDistances(aggregate[1], aggregate); - - map::iterator it; - set intersection; - for (it = a->begin(); it != a->end(); it++) - if (b->count((*it).first)) - intersection.insert((*it).first); - - // Continuing on for every envelope in the aggregate: - for (int i = 2; i < aggregate.size(); i++) - { - a = FindDistances(aggregate[i], aggregate); - set toRemove; - set::iterator sIt; - for (sIt = intersection.begin(); sIt != intersection.end(); sIt++) - if(a->count(*sIt) == 0) - toRemove.insert(*sIt); - for (sIt = toRemove.begin(); sIt != toRemove.end(); sIt++) - intersection.erase(*sIt); - } - - // Now finding the enclosing sphere sizes for all nodes in intersection - set::iterator sIt; - int bestSphereSize = 1000000; - for(sIt = intersection.begin(); sIt != intersection.end(); sIt++) - { - a = FindDistances(*sIt, aggregate); - int farthestInternal = 0; - for (it = a->begin(); it != a->end(); it++) - if (std::binary_search(aggregate.begin(), aggregate.end(), (*it).first)) - farthestInternal = std::max((*it).second, farthestInternal); - - int nodeCount = 0; - for (it = a->begin(); it != a->end(); it++) - if ((*it).second <= farthestInternal) - nodeCount++; - - if (bestSphereSize > nodeCount) - { - bestSphereSize = nodeCount; - } - } - - // Returning the ratio score - - return (double)aggregate.size()/bestSphereSize; - } - void MetricsContext::MakeConvex(vector& aggregate) { - makeConvexCalls++; - // Set to keep track of 'must have' external nodes - set neededExternal; - - // Keep track of all options for other external nodes - vector > > externalPossibilities; - - // Find paths - for (int rootIdx = 0; rootIdx < aggregate.size(); rootIdx++) - { - int startNode = aggregate[rootIdx]; - map *distances = FindDistances(startNode, aggregate); - - // Get paths for each distinct pair of nodes in the aggregate - for (int endIdx = rootIdx + 1; endIdx < aggregate.size(); endIdx++) - { - int endNode = aggregate[endIdx]; - - vector > *paths = GetShortestPaths(startNode, endNode, *distances); - vector > *externals = FindExternalsInPaths(aggregate, paths); - - // If path not satisfied add to externalPossibilities - if (!IsPathSatisfied(neededExternal, *externals)) - externalPossibilities.push_back(vector >(*externals)); - - // Memory cleanup - delete paths; - delete externals; - } - } - - // Simplify until no longer possible: - bool allDone = false; - while (!allDone) - { - allDone = true; - int counter = 0; - while (counter < externalPossibilities.size()) - { - if (IsPathSatisfied(neededExternal, externalPossibilities[counter])) - { - externalPossibilities.erase(externalPossibilities.begin() + counter); - allDone = false; - } - else - { - counter++; - } - } - } - - // If there were required externals found add them and recurse: - if (neededExternal.size() > 0) - { - - aggregate.insert(aggregate.end(), neededExternal.begin(), neededExternal.end()); - std::sort(aggregate.begin(), aggregate.end()); - MakeConvex(aggregate); - } - - // If there are external possibilities and no required we need brute force: - else if (externalPossibilities.size() > 0) - { - set* result = BruteForceMinimalNodes(externalPossibilities); - - - aggregate.insert(aggregate.end(), result->begin(), result->end()); - std::sort(aggregate.begin(), aggregate.end()); - delete result; - MakeConvex(aggregate); - } - - // Method falls through when the aggregate is convex - } - void MetricsContext::EnsureConvex(int aggId) { - if (aggId < aggregates.size() && aggId >= 0) - { - if (convexAggregates[aggId].size() < aggregates[aggId].size()) - { - convexAggregates[aggId].assign(aggregates[aggId].begin(), aggregates[aggId].end()); - MakeConvex(convexAggregates[aggId]); - } - } - } - vector* MetricsContext::FindCentroid(vector& aggregate) { - // To store the mass scores in: - map scores; - int currentNode = aggregate[0]; - - // Find score for first node - int bestScore = FindMassScore(currentNode, aggregate); - scores[currentNode] = bestScore; - - bool betterFound = true; - while(betterFound) - { - betterFound = false; - int start = graph->adjIndexes[currentNode]; - int end = graph->adjIndexes[currentNode + 1]; - for (int i = start; i < end && !betterFound; i++) - { - int neighbor = graph->adjacency[i]; - if (scores.count(neighbor) == 0) - scores[neighbor] = FindMassScore(neighbor, aggregate); - - if (scores[neighbor] < bestScore) - { - bestScore = scores[neighbor]; - currentNode = neighbor; - betterFound = true; - } - } - } - - // Find any adjacent nodes with equivalent score - vector *result = new vector(); - result->push_back(currentNode); - int start = graph->adjIndexes[currentNode]; - int end = graph->adjIndexes[currentNode + 1]; - for (int i = start; i < end; i++) - { - int neighbor = graph->adjacency[i]; - if (scores.count(neighbor) == 0) - scores[neighbor] = FindMassScore(neighbor, aggregate); - - if (scores[neighbor] == bestScore) - { - result->push_back(neighbor); - } - } - return result; - } - int MetricsContext::FindMassScore(int node, vector& aggregate) { - map *distances = FindDistances(node, aggregate); - int score = 0; - for (int i = 0; i < aggregate.size(); i++) - score += (*distances)[aggregate[i]]; - delete distances; - return score; - } - map* MetricsContext::FindDistances(int rootNode, vector& aggregate) { - distanceLookups++; - // Create queue and set distance of initial node - map *distances = new map(); - (*distances)[rootNode] = 0; - - // Queue for exploring out from root node - queue exploring; - exploring.push(rootNode); - - // Keeping track of whether the exploration should continue: - int internalsSeen = 0; // How many nodes in the aggregate have been marked - - // Start to explore the queue: - while (!exploring.empty()) - { - // Pull node off of queue - int node = exploring.front(); - exploring.pop(); - - // Get distance of current node: - int currentDistance = (*distances)[node]; - - // Examine the neighbors - for (int j = graph->adjIndexes[node]; j < graph->adjIndexes[node + 1]; j++) - { - int neighbor = graph->adjacency[j]; - if (distances->count(neighbor) == 0) - { - (*distances)[neighbor] = currentDistance + 1; - exploring.push(neighbor); - } - } - - // Checking if the current node is an internal: - if (std::binary_search(aggregate.begin(), aggregate.end(), node)) - { - internalsSeen++; - } - - // If we have seen all nodes in the aggregate stop the search - if (internalsSeen == aggregate.size()) - { - while(!exploring.empty()) - exploring.pop(); - break; - } - } - return distances; - } - vector >* MetricsContext::GetShortestPaths(int startId, int endId, map &distances) { - vector >* result = new vector >(1); - vector > &paths = *result; - - // Starting the trace back with current node: - paths[0].push_back(endId); - int activePath = 0; - while (activePath < paths.size()) - { - // Get the current end of the current path and its distance - int endNode = paths[activePath].back(); - - int distance = distances[endNode]; - bool branched = false; - - // Check the neighbors of the current endNode to continue - for (int nIt = graph->adjIndexes[endNode]; nIt < graph->adjIndexes[endNode + 1]; nIt++) - { - // Get neighbor and its distance - int neighbor = graph->adjacency[nIt]; - int neighborDist = distances[neighbor]; - - // If the neighbor is one closer than current it is on - // a shortest path - if (neighborDist == distance - 1) - { - // If there is a branch add a new path - if (branched) - { - // Add a new path starting with a copy of current - // and append current neighbor - paths.push_back(vector(paths[activePath].begin(),paths[activePath].end() - 1)); - paths.back().push_back(neighbor); - } - // Add node to current path and mark as branched - else - { - paths[activePath].push_back(neighbor); - branched = true; - } - } - } - - // Check if the active path is complete and move on - if (paths[activePath].back() == startId) - activePath++; - } - return result; - } - vector >* MetricsContext::FindExternalsInPaths(vector& aggregate, vector >* p) { - // Create the result vector and get a local reference to the paths - vector > *result = new vector >(); - vector > &paths = *p; - - // Go through each path and find external nodes to add - for (int p = 0; p < paths.size(); p++) - { - vector externals; - for (int pp = 0; pp < paths[p].size(); pp++) - { - int pathNode = paths[p][pp]; - - if (!std::binary_search(aggregate.begin(), aggregate.end(), pathNode)) - { - externals.push_back(pathNode); - } - } - - // If there were external nodes found add to the list - if (externals.size() > 0) - { - result->push_back(vector(externals)); - } - - // If a clear path was found return an empty path list - else - { - result->clear(); - return result; - } - } - - // Only keep meaningful elements - vector keepers(result->size(), true); - vector > &r = *result; - for(int a = 0; a < r.size(); a++) - { - // Compare with each following node: - for (int b = a + 1; b < r.size() && keepers[a]; b++) - { - // Make sure the b node is also not marked for deletion - if (keepers[b]) - { - // If the nodes are equal eliminate one - if (r[a].size() == r[b].size()) - { - // Assume that they are matched: - bool matched = true; - for (int i = 0; i < r[a].size(); i++) - { - bool good = false; - for (int j = 0; j < r[b].size(); j++) - { - if (r[a][i] == r[b][j]) - { - good = true; - break; - } - } - if (!good) - { - matched = false; - break; - } - } - if (matched) - { - keepers[b] = false; - } - } - - // Otherwise if the every element in the smaller is in the - // larger as well, eliminate the larger. - else - { - int small = a; - int big = b; - if (a > b) - { - small = b; - big = a; - } - - // Assume that they are matched: - bool matched = true; - for (int i = 0; i < r[small].size(); i++) - { - bool good = false; - for (int j = 0; j < r[big].size(); j++) - { - if (r[small][i] == r[big][j]) - { - good = true; - break; - } - } - if (!good) - { - matched = false; - break; - } - } - if (matched) - { - keepers[big] = false; - } - } - } - } - } - int toRemove = std::count(keepers.begin(), keepers.end(), false); - - if (toRemove == 0) - return result; - - // If there were some to remove do it: - vector > *trimmed = new vector >(); - trimmed->resize(r.size() - toRemove); - int insertNext = 0; - for (int i = 0; i < r.size(); i++) - if (keepers[i]) - (*trimmed)[insertNext++] = r[i]; - - delete result; - return trimmed; - } - bool MetricsContext::IsPathSatisfied(set& required, vector >& pathOptions) { - // If there are no options it is satisfied. - if (pathOptions.size() == 0) - return true; - - // If there is only one option those nodes are required so add and be done: - if (pathOptions.size() == 1) - { - required.insert(pathOptions[0].begin(), pathOptions[0].end()); - return true; - } - - // If any option is satisfied by nodes already required we are done: - for (int i = 0; i < pathOptions.size(); i++) - { - bool satisfied = true; - for (int j = 0; j < pathOptions[i].size(); j++) - if (required.count(pathOptions[i][j]) == 0) - satisfied = false; - if (satisfied) - { - return true; - } - } - - return false; - } - set* MetricsContext::BruteForceMinimalNodes(vector > >& pathOptions) { - set attempt; - for (int cIt = 0; cIt < pathOptions.size(); cIt++) - attempt.insert(pathOptions[cIt][0].begin(),pathOptions[cIt][0].end()); - - // Finding how many possibilites there are to check: - int possibilities = pathOptions[0].size(); - for (int i = 1; i < pathOptions.size(); i++) - possibilities *= pathOptions[i].size(); - if (possibilities > 100) - { - printf("%d options for brute forcing!\n", possibilities); - debugHelpers::printResult(&pathOptions, "Options"); - - } - - // Set best found to initial attempt size - int bestFound = attempt.size(); - set bestFoundSet = attempt; - - // Setting the guess - vector guess(pathOptions.size(), 0); - - // Trying all remaining combinations - int guessCount = 1; - while(IncrementGuessVector(guess, pathOptions)) - { - guessCount++; - - // Clearing attempt to known nodes - attempt.clear(); - - // Making choices according to guess vector - for (int cIt = 0; cIt < pathOptions.size(); cIt++) - { - attempt.insert(pathOptions[cIt][0].begin(),pathOptions[cIt][0].end()); - } - - // If this attempt has better result save it - if (attempt.size() < bestFound) - { - bestFound = attempt.size(); - bestFoundSet = attempt; - } - } - - set *result = new set(bestFoundSet); - return result; - } - bool MetricsContext::IncrementGuessVector(vector& guess, vector > >& externalOptions) { - bool incremented = false; - if (guess.size() == 0) - return false; - - // Increment the guess vector - int position = guess.size() - 1; - while (true) - { - // If the current position has a higher option increment and stop: - if (guess[position] < externalOptions[position].size() - 1) - { - guess[position]++; - incremented = true; - break; - } - // If the current position can go no higher: - else - { - // If there is no preceding position all options have been explored. - if (position == 0) - break; - - // Otherwise reset this position to zeros and move up to next position. - else - { - guess[position] = 0; - position--; - } - } - } - return incremented; - } - - // Setup Helpers - void MetricsContext::Initialize() { - GetAggregates(); - convexAggregates.resize(aggregates.size()); - } - void MetricsContext::GetAggregates() { - // Iterate through the aggregation and populate aggregates - for (int i = 0; i < aggregation.size(); i++) - { - int aggId = aggregation[i]; - - // Make sure an entry for this aggregate exists - if (aggregates.size() <= aggId) - { - aggregates.resize(aggId + 1); - } - - // Push the current node onto the list for its aggregate. - aggregates[aggId].push_back(i); - } - } - } -} diff --git a/src/core/aggmis/cuda/AggMIS_Types.cu b/src/core/aggmis/cuda/AggMIS_Types.cu deleted file mode 100644 index 7b02fb3..0000000 --- a/src/core/aggmis/cuda/AggMIS_Types.cu +++ /dev/null @@ -1,432 +0,0 @@ -/* - * File: AggMIS_Types.cu - * Author: T. James Lewis - * - * Created on April 15, 2013, 2:18 PM - */ -#include "AggMIS_Types.h" - -namespace AggMIS { - bool CheckCudaError(cudaError_t code, const char* file, int line) { - if (code != cudaSuccess) { - std::cout << "\n***************** CUDA Error detected ***************\n"; - std::cout << "Error: " << cudaGetErrorString(code) << "\n"; - std::cout << "In file " << file << " line " << line << "\n"; - std::cout << "\n*****************************************************\n"; - } - code = cudaGetLastError(); - if (code != cudaSuccess) { - std::cout << "\n*************** Past CUDA Error detected ************\n"; - std::cout << "Error: " << cudaGetErrorString(code) << "\n"; - std::cout << "In file " << file << " line " << line << "\n"; - std::cout << "\n*****************************************************\n"; - } - return false; - } - namespace Types { - // My timer implementation - JTimer::JTimer() { - cudaEventCreate(&startTimeCuda); - cudaEventCreate(&endTimeCuda); - started = false; - stopped = false; - startTimeHost.tv_sec = startTimeHost.tv_nsec = 0; - endTimeHost.tv_sec = endTimeHost.tv_nsec = 0; - } - JTimer::~JTimer() {} - void JTimer::start() { - cudaEventRecord(startTimeCuda, 0); - clock_gettime(CLOCK_REALTIME, &startTimeHost); - started = true; - stopped = false; - } - void JTimer::stop() { - if (started && !stopped) { - cudaEventRecord(endTimeCuda, 0); - cudaEventSynchronize(endTimeCuda); - clock_gettime(CLOCK_REALTIME, &endTimeHost); - stopped = true; - } - } - double JTimer::getElapsedTimeInSec(bool host) { - if (!started || !stopped) { - printf("Error: elapsed time requested when not valid.\n"); - return -1.0; - } - if (host) { - if (endTimeHost.tv_nsec < startTimeHost.tv_nsec) { - endTimeHost.tv_nsec += 1000000000; - endTimeHost.tv_sec -= 1; - } - long timeInMicrosec = ((endTimeHost.tv_sec - startTimeHost.tv_sec) * 1000000) - + ((endTimeHost.tv_nsec - startTimeHost.tv_nsec) / 1000); - return (double)(timeInMicrosec) / 1000000.0; - } - else { - cudaEventElapsedTime(&elapsedCudaTime, startTimeCuda, endTimeCuda); - return (double) elapsedCudaTime / 1000.0; - } - } - double JTimer::getElapsedTimeInMilliSec(bool host) { - if (host) { - // Check if we need to carry some nanoseconds - if (endTimeHost.tv_nsec < startTimeHost.tv_nsec) { - endTimeHost.tv_nsec += 1000000000; - endTimeHost.tv_sec -= 1; - } - long timeInMicrosec = ((endTimeHost.tv_sec - startTimeHost.tv_sec) * 1000) - + ((endTimeHost.tv_nsec - startTimeHost.tv_nsec) / 1000); - return (double)(timeInMicrosec) / 1000.0; - } - else { - cudaEventElapsedTime(&elapsedCudaTime, startTimeCuda, endTimeCuda); - return (double) elapsedCudaTime; - } - } - - // Graph_d members - Graph_d::Graph_d(IntVector_d& indices, - IntVector_d& adjacency) { - this->indices = new IntVector_d(indices); - this->adjacency = new IntVector_d(adjacency); - willClean = true; - } - Graph_d::Graph_d(IntVector_h& indices, - IntVector_h& adjacency) { - this->indices = new IntVector_d(indices); - this->adjacency = new IntVector_d(adjacency); - willClean = true; - } - Graph_d::Graph_d(IntVector_d* indices, - IntVector_d* adjacency) { - this->indices = indices; - this->adjacency = adjacency; - willClean = false; - } - Graph_d::Graph_d(Graph_h& graph) { - indices = new IntVector_d(*(graph.indices)); - adjacency = new IntVector_d(*(graph.adjacency)); - willClean = true; - } - Graph_d::Graph_d() { - indices = new IntVector_d(); - adjacency = new IntVector_d(); - willClean = true; - } - Graph_d::~Graph_d() { - if (willClean) - { - indices->clear(); - adjacency->clear(); - delete indices; - delete adjacency; - } - } - int Graph_d::Size() { - return indices->size() - 1; - } - DGraph Graph_d::GetD() { - return DGraph(Size(), indStart(), adjStart()); - } - int* Graph_d::indStart() { - return thrust::raw_pointer_cast(indices->data()); - } - int* Graph_d::adjStart() { - return thrust::raw_pointer_cast(adjacency->data()); - } - - // Graph_h members - Graph_h::Graph_h(IntVector_d& indices, - IntVector_d& adjacency) { - this->indices = new IntVector_h(indices); - this->adjacency = new IntVector_h(adjacency); - willClean = true; - } - Graph_h::Graph_h(IntVector_h& indices, - IntVector_h& adjacency) { - this->indices = new IntVector_h(indices); - this->adjacency = new IntVector_h(adjacency); - willClean = true; - } - Graph_h::Graph_h(IntVector_h* indices, - IntVector_h* adjacency) { - this->indices = indices; - this->adjacency = adjacency; - willClean = false; - } - Graph_h::Graph_h(Graph_d& graph) { - indices = new IntVector_h(*(graph.indices)); - adjacency = new IntVector_h(*(graph.adjacency)); - willClean = true; - } - Graph_h::Graph_h() { - indices = new IntVector_h(); - adjacency = new IntVector_h(); - willClean = true; - } - Graph_h::~Graph_h() { - if (willClean) - { - indices->resize(0); - adjacency->resize(0); - delete indices; - delete adjacency; - } - } - int Graph_h::Size() { - return indices->size() - 1; - } - int* Graph_h::nStart(int node) { - return &((*adjacency)[(*indices)[node]]); - } - int* Graph_h::nEnd(int node) { - return &((*adjacency)[(*indices)[node + 1]]); - } - - // Functions - int* StartOf(IntVector_d &target) { - return thrust::raw_pointer_cast(target.data()); - } - int* StartOf(IntVector_d *target) { - return thrust::raw_pointer_cast(target->data()); - } - - namespace Compare { - bool AreEqual(IntVector_h& a, - IntVector_h& b, - bool verbose) { - bool good = true; - if (a.size() != b.size()) - { - if (verbose) - printf("Vectors to compare differ in size: a.size()=%d b.size=%d\n", a.size(), b.size()); - return false; - } - - for (int i = 0; i < a.size(); i++) - if (a[i] != b[i]) - { - if (verbose) - printf("Difference found: a[%d]=%d b[%d]=%d\n", i, a[i], i, b[i]); - good = false; - } - return good; - } - bool AreEqual(IntVector_d& a, - IntVector_d& b, - bool verbose) { - IntVector_h tempA(a); - IntVector_h tempB(b); - bool result = AreEqual(tempA, tempB, verbose); - tempA.clear(); - tempB.clear(); - return result; - } - bool AreEqual(IntVector_h& a, - IntVector_d& b, - bool verbose) { - IntVector_h temp(b); - bool result = AreEqual(a, temp, verbose); - temp.clear(); - return result; - } - bool AreEqual(IntVector_d& a, - IntVector_h& b, - bool verbose) { - return AreEqual(b, a, verbose); - } - bool AreEqual(vector > &a, - vector > &b, - bool verbose) { - // Check that main containers have matching sizes - if (a.size() != b.size()) { - if (verbose) - printf("Sizes of base vectors do not match! a=%d b=%d\n", - a.size(), b.size()); - return false; - } - // Check that sizes of nested containers match - for (int i = 0; i < a.size(); i++) { - if (a[i].size() != b[i].size()) { - if (verbose) { - printf("Sizes of secondary vectors %d do not match!\n", i); - printf("a[%d].size()=%d b[%d].size()=%d\n", - i, a[i].size(), i, b[i].size()); - stringstream ss; - ss << "Contents of A[" << i << "]"; - Display::Print(a[i], ss.str()); - ss.str("Contents of B["); - ss << i << "]"; - Display::Print(b[i], ss.str()); - } - return false; - } - } - // Check that all entries are equal - for (int i = 0; i < a.size(); i++) { - for (int j = 0; j < a[i].size(); j++) { - if (a[i][j] != b[i][j]) { - if (verbose) { - printf("Element[%d][%d] does not match!\n", i, j); - stringstream ss; - ss << "Contents of A[" << i << "]"; - Display::Print(a[i], ss.str()); - ss.str("Contents of B["); - ss << i << "]"; - Display::Print(b[i], ss.str()); - } - return false; - } - } - } - return true; - } - bool AreEqual(Graph_h& a, - Graph_h& b, - bool verbose) { - bool indicesMatch = AreEqual(*(a.indices), - *(b.indices), - verbose); - bool adjacencyMatch = AreEqual(*(a.adjacency), - *(b.adjacency), - verbose); - if (!indicesMatch && verbose) - printf("Indices of graphs differ!\n"); - if (!adjacencyMatch && verbose) - printf("Adjacency lists of graphs differ!\n"); - return indicesMatch && adjacencyMatch; - } - bool AreEqual(Graph_d& a, - Graph_d& b, - bool verbose) { - bool indicesMatch = AreEqual(*(a.indices), - *(b.indices), - verbose); - bool adjacencyMatch = AreEqual(*(a.adjacency), - *(b.adjacency), - verbose); - if (!indicesMatch && verbose) - printf("Indices of graphs differ!\n"); - if (!adjacencyMatch && verbose) - printf("Adjacency lists of graphs differ!\n"); - return indicesMatch && adjacencyMatch; - } - bool AreEqual(Graph_h& a, - Graph_d& b, - bool verbose) { - bool indicesMatch = AreEqual(*(a.indices), - *(b.indices), - verbose); - bool adjacencyMatch = AreEqual(*(a.adjacency), - *(b.adjacency), - verbose); - if (!indicesMatch && verbose) - printf("Indices of graphs differ!\n"); - if (!adjacencyMatch && verbose) - printf("Adjacency lists of graphs differ!\n"); - return indicesMatch && adjacencyMatch; - } - bool AreEqual(Graph_d& a, - Graph_h& b, - bool verbose) { - bool indicesMatch = AreEqual(*(a.indices), - *(b.indices), - verbose); - bool adjacencyMatch = AreEqual(*(a.adjacency), - *(b.adjacency), - verbose); - if (!indicesMatch && verbose) - printf("Indices of graphs differ!\n"); - if (!adjacencyMatch && verbose) - printf("Adjacency lists of graphs differ!\n"); - return indicesMatch && adjacencyMatch; - } - } - namespace Display { - void Print(IntVector_h& toPrint, - int start, - int end, - string message) { - printf("%s:\n", message.c_str()); - printf("\n %8d: ", 0); - for (int i = start; i < end; i++) - { - if ((i-start) % 10 == 0 && (i-start) > 0) - printf("\n %8d: ", i); - - int value = toPrint[i]; - printf(" %8d", value); - } - printf("\n"); - } - void Print(IntVector_d& toPrint, - int start, - int end, - string message) { - IntVector_h temp(toPrint); - Print(temp, start, end, message); - temp.clear(); - } - void Print(IntVector_d& toPrint, - string message) { - IntVector_h temp(toPrint); - Print(temp, 0, temp.size(), message); - temp.clear(); - } - void Print(IntVector_h& toPrint, - string message) { - Print(toPrint, 0, toPrint.size(), message); - } - void Print(vector > >& toPrint, - string message) { - // Print out general info: - printf("Triple vector %s has %d entries:\n", message.c_str(), toPrint.size()); - - for (int i = 0; i < toPrint.size(); i++) - { - cout << message << "[" << i << "]: "; - for (int z = 0; z < toPrint[i].size(); z++) - { - cout << "("; - for (int zz = 0; zz < toPrint[i][z].size(); zz++) - { - cout << toPrint[i][z][zz]; - if (zz < toPrint[i][z].size() -1) - cout << " "; - } - cout << ") "; - } - cout << "\n"; - } - cout << "\n"; - } - void Print(vector >& toPrint, - string message) { - printf("%s:\n", message.c_str()); - for (int j = 0; j < toPrint.size(); j++) { - printf("\n %4d: ", j); - for (int i = 0; i < toPrint[j].size(); i++) - { - if (i % 10 == 0 && i > 0) - printf("\n %4d: ", j); - - int value = toPrint[j][i]; - printf(" %4d", value); - } - } - printf("\n"); - } - void Print(vector& toPrint, - int start, - int end, - string message) { - IntVector_h temp(toPrint.begin(), toPrint.end()); - Print(temp, start, end, message); - } - void Print(vector& toPrint, - string message) { - Print(toPrint, 0, toPrint.size(), message); - } - } - } -} diff --git a/src/core/cuda/setup_solver.cu b/src/core/cuda/setup_solver.cu index b1074a7..b47558f 100644 --- a/src/core/cuda/setup_solver.cu +++ b/src/core/cuda/setup_solver.cu @@ -512,3 +512,36 @@ int writeMatlabArray(const std::string &filename, const Vector_h_CG &array) { return 0; } + +void writeVTK(TetMesh* m_meshPtr, std::vector values, std::string fname) +{ + int nv = m_meshPtr->vertices.size(); + int nt = m_meshPtr->tets.size(); + FILE* vtkfile; + vtkfile = fopen((fname + ".vtk").c_str(), "w+"); + fprintf(vtkfile, "# vtk DataFile Version 3.0\nvtk output\nASCII\nDATASET UNSTRUCTURED_GRID\n"); + fprintf(vtkfile, "POINTS %d float\n", nv); + for (int i = 0; i < nv; i++) + { + fprintf(vtkfile, "%.12f %.12f %.12f\n", m_meshPtr->vertices[i][0], + m_meshPtr->vertices[i][1], m_meshPtr->vertices[i][2]); + } + fprintf(vtkfile, "CELLS %d %d\n", nt, nt * 5); + for (int i = 0; i < nt; i++) + { + fprintf(vtkfile, "4 %d %d %d %d\n", m_meshPtr->tets[i][0], + m_meshPtr->tets[i][1], m_meshPtr->tets[i][2], m_meshPtr->tets[i][3]); + } + + fprintf(vtkfile, "CELL_TYPES %d\n", nt); + for (int i = 0; i < nt; i++) + { + fprintf(vtkfile, "10\n"); + } + fprintf(vtkfile, "POINT_DATA %d\nSCALARS traveltime float 1\nLOOKUP_TABLE default\n", + nv, values.size()); + for (int i = 0; i < values.size(); i++) { + fprintf(vtkfile, "%.12f\n ", values[i]); + } + fclose(vtkfile); +} \ No newline at end of file diff --git a/src/core/include/FEM.h b/src/core/include/FEM.h new file mode 100644 index 0000000..2865acb --- /dev/null +++ b/src/core/include/FEM.h @@ -0,0 +1,147 @@ +#ifndef __FEM_H__ +#define __FEM_H__ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace FEM { + /** The class that represents all of the available options for FEM */ + class FEM { + public: + FEM(std::string fname = "../src/test/test_data/sphere334", + bool verbose = false) : + verbose_(verbose), + filename_(fname), + seedPointList_(std::vector(1, 0)), + maxBlocks_(1000), + maxVertsPerBlock_(64), + isStructured_(false), + squareLength_(16), + squareWidth_(16), + squareDepth_(16), + squareBlockLength_(1), + squareBlockWidth_(1), + squareBlockDepth_(1), + maxIterations_(1000) + + {} + //3D data + bool verbose_; + std::string filename_; + std::vector seedPointList_; + int maxBlocks_; + int maxVertsPerBlock_; + bool isStructured_; + int squareLength_, squareWidth_, squareDepth_; + int squareBlockLength_, squareBlockWidth_, squareBlockDepth_; + int maxIterations_; + }; + + //The static pointer to the mesh + static TetMesh * mesh_ = NULL; + //the answer vector + static std::vector < std::vector > iteration_values_; + //accessor functions to the results. + std::vector < float >& getFinalResult() { + return iteration_values_.at(iteration_values_.size() - 1); + } + std::vector < float >& getResultAtIteration(size_t i) { + return iteration_values_.at(i); + } + size_t numIterations() { return iteration_values_.size(); } + void writeVTK() { + //TODO + } + + /** + * Creates the mesh, partitions the mesh, and runs the algorithm. + * + * @data The set of options for the Eikonal algorithm. + * The defaults are used if nothing is provided. + */ + void solveFEM(FEM data = FEM()) { + + clock_t starttime, endtime; + starttime = clock (); + mesh_ = TetMesh::read( + (data.filename_ + ".node").c_str(), + (data.filename_ + ".ele").c_str(), true, data.verbose_); + mesh_->rescale(4.0); + mesh_->need_neighbors(data.verbose_); + mesh_->need_meshquality(data.verbose_); + //TODO + endtime = clock(); + double duration = (double)(endtime - starttime) * 1000/ CLOCKS_PER_SEC; + + if (data.verbose_) + printf("Computing time : %.10lf ms\n",duration); + } + + /** + * This function uses the provided analytical solutions to + * visualize the algorithm's error after each iteration. + * + * @param solution The vector of expected solutions. + */ + void printErrorGraph(std::vector solution) { + + // now calculate the RMS error for each iteration + std::vector rmsError; + rmsError.resize(numIterations()); + for (size_t i = 0; i < numIterations(); i++) { + float sum = 0.f; + std::vector result = getResultAtIteration(i); + for (size_t j = 0; j < solution.size(); j++) { + float err = std::abs(solution[j] - result[j]); + sum += err * err; + } + rmsError[i] = std::sqrt(sum / static_cast(solution.size())); + } + //determine the log range + float max_err = rmsError[0]; + float min_err = rmsError[rmsError.size() - 1]; + int max_log = -10, min_log = 10; + while (std::pow(static_cast(10),max_log) < max_err) max_log++; + while (std::pow(static_cast(10),min_log) > min_err) min_log--; + // print the error graph + printf("\n\nlog(Err)|\n"); + bool printTick = true; + for(int i = max_log ; i >= min_log; i--) { + if (printTick) { + printf(" 10^%2d|",i); + } else { + printf(" |"); + } + for (size_t j = 0; j < numIterations(); j++) { + if (rmsError[j] > std::pow(static_cast(10),i) && + rmsError[j] < std::pow(static_cast(10),i+1)) + printf("*"); + else + printf(" "); + } + printf("\n"); + printTick = !printTick; + } + printf("--------|------------------------------------------"); + printf(" Converged to: %.4f\n",rmsError[rmsError.size() - 1]); + printf(" |1 5 10 15 20 25 30 35\n"); + printf(" Iteration\n"); + } +} + +#endif diff --git a/src/core/include/my_timer.h b/src/core/include/my_timer.h index 2d05463..2aba38f 100644 --- a/src/core/include/my_timer.h +++ b/src/core/include/my_timer.h @@ -26,6 +26,7 @@ int inline clock_gettime(int clk_id, struct timespec *t){ double inline CLOCK() { #ifdef WIN32 #include +#include SYSTEMTIME st; GetSystemTime(&st); return ((st.wDay * 24. + st.wHour) * 60. + st.wMinute) * 60. + st.wSecond + st.wMilliseconds / 1000.; diff --git a/src/core/include/setup_solver.h b/src/core/include/setup_solver.h index e1b53a0..d3cafd6 100644 --- a/src/core/include/setup_solver.h +++ b/src/core/include/setup_solver.h @@ -53,4 +53,7 @@ int writeMatlabArray(const std::string &filename, const Vector_h_CG &array); void checkMatrixForValidContents(Matrix_ell_h* A_h, const bool verbose); + +void writeVTK(TetMesh* m_meshPtr, std::vector values, std::string fname); + #endif diff --git a/src/examples/example1.cu b/src/examples/example1.cu index 54094b0..4798897 100644 --- a/src/examples/example1.cu +++ b/src/examples/example1.cu @@ -82,6 +82,8 @@ int main(int argc, char** argv) getMatrixFromMesh(cfg, tetmeshPtr, &A_h, verbose); //intialize the b matrix to ones for now. TODO @DEBUG Vector_h_CG b_h(A_h.num_rows, 1.0); + for (size_t i = 0; i < b_h.size() / 2; i++) + b_h[i] = 0.; //The answer vector. Vector_h_CG x_h(A_h.num_rows, 0.0); //intial X vector //************************ DEBUG : creating identity matrix for stiffness properties for now. @@ -108,5 +110,15 @@ int main(int argc, char** argv) if (writeMatlabArray("output.mat", x_h)) { std::cerr << "failed to write matlab file." << std::endl; } + //write the VTK + std::vector vals; + for (size_t i = 0; i < x_h.size(); i++){ + vals.push_back(x_h[i]); + } + auto pos = filename.find_last_of("/"); + if (pos == std::string::npos) + pos = filename.find_last_of("\\"); + filename = filename.substr(pos + 1, filename.size() - 1); + writeVTK(tetmeshPtr, vals, filename); return 0; }