Brandon Potter
1600 строки
69 KiB
C++
1600 строки
69 KiB
C++
/********************************************************************************
|
|
* Copyright (c) 2024 Advanced Micro Devices, Inc. All rights reserved.
|
|
*
|
|
* Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
* of this software and associated documentation files (the "Software"), to deal
|
|
* in the Software without restriction, including without limitation the rights
|
|
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
* copies of the Software, and to permit persons to whom the Software is
|
|
* furnished to do so, subject to the following conditions:
|
|
*
|
|
* The above copyright notice and this permission notice shall be included in
|
|
* all copies or substantial portions of the Software.
|
|
*
|
|
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
|
|
* THE SOFTWARE.
|
|
********************************************************************************/
|
|
#ifndef SpTS_H
|
|
#define SpTS_H
|
|
#define TEST_NUM 9999999999ULL
|
|
|
|
#include "InputFlags.h"
|
|
#include "SparseMatrix.h"
|
|
|
|
#include "MatrixMarketReader.h"
|
|
#include <vector>
|
|
#include <float.h>
|
|
#include <unordered_map>
|
|
#include <cmath>
|
|
#include <chrono>
|
|
#include <thread>
|
|
|
|
#include <time.h>
|
|
#include <sys/time.h>
|
|
|
|
#include <unistd.h>
|
|
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
#include "roc_shmem.hpp"
|
|
#include "mpi.h"
|
|
#endif
|
|
|
|
#ifdef USE_HIP
|
|
#include "spts_kernel.h"
|
|
#endif
|
|
|
|
#ifdef DBL_DECIMAL_DIG
|
|
#define OP_DBL_Digs (DBL_DECIMAL_DIG)
|
|
#else
|
|
#ifdef DECIMAL_DIG
|
|
#define OP_DBL_Digs (DECIMAL_DIG)
|
|
#else
|
|
#define OP_DBL_Digs (DBL_DIG + 3)
|
|
#endif
|
|
#endif
|
|
|
|
#ifdef FLT_DECIMAL_DIG
|
|
#define OP_FLT_Digs (FLT_DECIMAL_DIG)
|
|
#else
|
|
#ifdef DECIMAL_DIG
|
|
#define OP_FLT_Digs (DECIMAL_DIG)
|
|
#else
|
|
#define OP_FLT_Digs (FLT_DIG + 3)
|
|
#endif
|
|
#endif
|
|
|
|
using namespace rocshmem;
|
|
|
|
template<typename FloatType>
|
|
class SparseTriangularSolve :
|
|
public InputFlags, public SparseMatrix<FloatType>
|
|
{
|
|
FloatType *x;
|
|
FloatType *y;
|
|
FloatType *y_zero;
|
|
FloatType *yref;
|
|
std::vector<uint64_t> rowBlocks;
|
|
|
|
memPointer xDev;
|
|
memPointer yDev;
|
|
memPointer completedRowsDev;
|
|
memPointer rowBlocksDev;
|
|
memPointer doneArrayDev;
|
|
memPointer shadowDoneArrayDev;
|
|
memPointer remoteInProgressArrayDev;
|
|
memPointer reqUpdateArrayDev;
|
|
memPointer numRowsAtLevelDev;
|
|
memPointer maxDepthDev;
|
|
memPointer rowMapDev;
|
|
memPointer totalSpinDev;
|
|
memPointer oneBufDev;
|
|
|
|
int nNZ;
|
|
int nRows;
|
|
int nCols;
|
|
int numBlocks;
|
|
/*
|
|
#ifdef USE_ROC_SHMEM
|
|
roc_shmem_t* handle;
|
|
#endif
|
|
*/
|
|
std::unordered_map<int, FloatType> *observed_errors;
|
|
int *errors_seen;
|
|
|
|
public:
|
|
|
|
SparseTriangularSolve() : nNZ(0), nRows(0), nCols(0), numBlocks(0)
|
|
{
|
|
x = NULL; y = NULL; y_zero = NULL, yref = NULL, observed_errors = NULL, errors_seen = NULL;
|
|
xDev = yDev = completedRowsDev = remoteInProgressArrayDev = rowBlocksDev = doneArrayDev = shadowDoneArrayDev = numRowsAtLevelDev = maxDepthDev = rowMapDev = totalSpinDev = oneBufDev = 0;
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
int roc_shmem_queues = (2560 / WF_PER_WG);
|
|
if (2560 % WF_PER_WG)
|
|
roc_shmem_queues++;
|
|
printf("roc_shmem_queues %d WF_PER_WG %d \n",roc_shmem_queues, WF_PER_WG);
|
|
roc_shmem_init(roc_shmem_queues);
|
|
|
|
this->Set_total_pes(roc_shmem_n_pes());
|
|
this->Set_this_pe(roc_shmem_my_pe());
|
|
#else
|
|
this->Set_total_pes(1);
|
|
this->Set_this_pe(0);
|
|
#endif
|
|
}
|
|
|
|
void AddDerivedInputFlags();
|
|
void AllocateVectors(MatrixMarketReader<FloatType> &mm_reader);
|
|
void CSRSpTSCPU(FloatType alpha);
|
|
bool CSRCheckCPU(FloatType alpha);
|
|
|
|
float CSRSpTSGPU(uint64_t &ns_per_iter, uint64_t &ns_per_analysis_iter, uint64_t &ns_per_syncfree_iter, uint64_t &ns_per_levelset_iter, uint64_t &ns_per_levelsync_iter, FloatType alpha);
|
|
|
|
int VerifyResults(int);
|
|
int NonDeterministicErrors();
|
|
int MaxErrors();
|
|
int ComputeRowBlocks(std::vector<uint64_t> &, int *, int);
|
|
|
|
~SparseTriangularSolve()
|
|
{
|
|
if (x != NULL)
|
|
delete[] x;
|
|
if (y != NULL)
|
|
delete[] y;
|
|
if (y_zero != NULL)
|
|
delete[] y_zero;
|
|
if (yref != NULL)
|
|
delete[] yref;
|
|
if (errors_seen != NULL)
|
|
delete[] errors_seen;
|
|
|
|
if (xDev != 0)
|
|
this->GPU->FreeMem(xDev);
|
|
if (rowBlocksDev != 0)
|
|
this->GPU->FreeMem(rowBlocksDev);
|
|
if (completedRowsDev != 0)
|
|
this->GPU->FreeMem(completedRowsDev);
|
|
if (numRowsAtLevelDev != 0)
|
|
this->GPU->FreeMem(numRowsAtLevelDev);
|
|
if (maxDepthDev != 0)
|
|
this->GPU->FreeMem(maxDepthDev);
|
|
if (rowMapDev != 0)
|
|
this->GPU->FreeMem(rowMapDev);
|
|
if (totalSpinDev != 0)
|
|
this->GPU->FreeMem(totalSpinDev);
|
|
if (oneBufDev != 0)
|
|
this->GPU->FreeMem(oneBufDev);
|
|
if (remoteInProgressArrayDev != 0)
|
|
this->GPU->FreeMem(remoteInProgressArrayDev);
|
|
|
|
#ifndef USE_ROC_SHMEM
|
|
if (yDev != 0)
|
|
this->GPU->FreeMem(yDev);
|
|
if (doneArrayDev != 0)
|
|
this->GPU->FreeMem(doneArrayDev);
|
|
if (reqUpdateArrayDev != 0)
|
|
this->GPU->FreeMem(reqUpdateArrayDev);
|
|
if (shadowDoneArrayDev != 0)
|
|
this->GPU->FreeMem(shadowDoneArrayDev);
|
|
#else
|
|
if (yDev != 0)
|
|
roc_shmem_free(yDev);
|
|
if (doneArrayDev != 0)
|
|
roc_shmem_free(doneArrayDev);
|
|
if (reqUpdateArrayDev != 0)
|
|
roc_shmem_free(reqUpdateArrayDev);
|
|
if (shadowDoneArrayDev != 0)
|
|
roc_shmem_free(shadowDoneArrayDev);
|
|
roc_shmem_finalize();
|
|
#endif
|
|
}
|
|
};
|
|
|
|
template<typename FloatType>
|
|
void SparseTriangularSolve<FloatType>::AddDerivedInputFlags()
|
|
{
|
|
AddInputFlag("filename", 'f', "", "Matrix-Market File", "string");
|
|
AddInputFlag("iterations", 'i', "10", "Number of SpTS Iterations (Default=10)", "int");
|
|
AddInputFlag("exp_zeroes", 'z', "false", "Include Explicit Zeroes in Matrix-Market File (Default=false)", "bool");
|
|
AddInputFlag("device", 'd', "0", "Choose the GPU to Execute SpTS (Default=0)", "int");
|
|
AddInputFlag("alpha", 'A', "1.0", "A*y=alpha*x. Known vector 'x' is multiplied by scalar alpha befoer solving for vector 'y'. (Default=1.0)", "float");
|
|
AddInputFlag("non_symmetric", 'n', "false", "Force the program to work on non-symmetric matrices. This will ignore the upper triangular entirely. (Default=false)", "bool");
|
|
AddInputFlag("levelsync_size", 'l', "0", "Number of rows to launch in a level-sync kernel invocation (Default = auto-tune)", "int");
|
|
AddInputFlag("verify", 'v', "false", "Verify results", "bool");
|
|
AddInputFlag("rocshmem_algorithm", 'a', "0", "ROC_SHMEM algorithm type", "int");
|
|
AddInputFlag("block_size", 'b', "32768", "Use get-based algorithm for ROC_SHMEM", "int");
|
|
AddInputFlag("put_block_size", 'p', "1024", "Block size for puts", "int");
|
|
AddInputFlag("get_backoff_factor", 'g', "128", "Backoff factor for gets", "int");
|
|
}
|
|
|
|
template<typename FloatType>
|
|
void SparseTriangularSolve<FloatType>::AllocateVectors(
|
|
MatrixMarketReader<FloatType> &mm_reader)
|
|
{
|
|
nRows = mm_reader.GetNumRows();
|
|
nCols = mm_reader.GetNumCols();
|
|
nNZ = mm_reader.GetNumNonZeroes();
|
|
|
|
x = new FloatType[nCols];
|
|
y = new FloatType[nRows];
|
|
y_zero = new FloatType[nRows];
|
|
yref = new FloatType[nRows];
|
|
observed_errors = new std::unordered_map<int, FloatType>[InputFlags::GetValueInt("iterations")];
|
|
|
|
for(int i = 0; i < nRows; i++)
|
|
{
|
|
y[i] = (FloatType)0.0;
|
|
y_zero[i] = (FloatType)0.0;
|
|
yref[i] = (FloatType)0.0;
|
|
}
|
|
|
|
for(int i = 0; i < nCols; i++)
|
|
{
|
|
//x[i] = (FloatType)rand() / (FloatType)RAND_MAX;
|
|
x[i] = 2.;
|
|
}
|
|
|
|
xDev = this->GPU->AllocateMem("xDev", nCols*sizeof(FloatType), GPU_MEM_READ_ONLY, NULL);
|
|
#ifndef USE_ROC_SHMEM
|
|
yDev = this->GPU->AllocateMem("yDev", nRows*sizeof(FloatType), GPU_MEM_READ_WRITE, NULL);
|
|
#else
|
|
yDev = (memPointer) roc_shmem_malloc(nRows*sizeof(FloatType));
|
|
#endif
|
|
}
|
|
|
|
template<typename FloatType>
|
|
void SparseTriangularSolve<FloatType>::CSRSpTSCPU(FloatType alpha)
|
|
{
|
|
FloatType *NZvalues = SparseMatrix<FloatType>::GetVals();
|
|
int *Cols = SparseMatrix<FloatType>::GetCols();
|
|
int *rowptrs = SparseMatrix<FloatType>::GetRowPtrs();
|
|
double internal_alpha = alpha;
|
|
|
|
uint64_t local_nnz = 0;
|
|
uint64_t remote_nnz = 0;
|
|
uint64_t rows_with_nonlocal = 0;
|
|
|
|
for(int i = 0; i < nRows; i++)
|
|
{
|
|
bool row_has_nonlocal = false;
|
|
double diagonal = 0.;
|
|
double temp = 0.;
|
|
int diag_j = -1;
|
|
for(int j = rowptrs[i]; j < rowptrs[i+1]; j++)
|
|
{
|
|
int ci = Cols[j];
|
|
int row_pe = (i / SPTS_BLOCK_SIZE) % this->Get_total_pes();
|
|
int col_pe = (ci / SPTS_BLOCK_SIZE) % this->Get_total_pes();
|
|
|
|
int assigned_pe = (i / SPTS_BLOCK_SIZE) % this->Get_total_pes();
|
|
if (assigned_pe == this->Get_this_pe()) {
|
|
if (row_pe == col_pe) {
|
|
local_nnz++;
|
|
} else {
|
|
row_has_nonlocal = true;
|
|
remote_nnz++;
|
|
}
|
|
}
|
|
|
|
// Skip adding in the diagonal. We need to solve for that.
|
|
if (ci != i)
|
|
{
|
|
if (i == TEST_NUM)
|
|
fprintf(stderr, "NZvalues[%d](%lf) * yref[%d](%lf)\n", j, NZvalues[j], ci, yref[ci]);
|
|
temp += NZvalues[j] * yref[ci];
|
|
}
|
|
else
|
|
{
|
|
if (i==TEST_NUM)
|
|
fprintf(stderr, "\t\tDIAG = %lf\n", NZvalues[j]);
|
|
diagonal = NZvalues[j];
|
|
diag_j = j;
|
|
}
|
|
}
|
|
if (row_has_nonlocal) rows_with_nonlocal++;
|
|
if (diag_j == -1)
|
|
{
|
|
fflush(stdout);
|
|
printf("\nERROR in SpTS CPU\n");
|
|
printf("No diagonal found in row %d\n", i);
|
|
}
|
|
// y = (x-sum_of_vals_from_A) / diag
|
|
double alpha_x = internal_alpha * (double)x[i];
|
|
if (i == TEST_NUM)
|
|
{
|
|
char buf[128];
|
|
char buf2[128];
|
|
char buf3[128];
|
|
snprintf(buf, sizeof(buf), "%.20f", alpha_x);
|
|
snprintf(buf2, sizeof(buf2), "%.20f", internal_alpha);
|
|
fprintf(stderr, "alpha_x: %s (%s * %lf)\n", buf, buf2, x[i]);
|
|
snprintf(buf3, sizeof(buf3), "%.20f", temp);
|
|
fprintf(stderr, "temp: %s\n", buf3);
|
|
}
|
|
yref[i] = (FloatType)((alpha_x - temp)/diagonal);
|
|
if (i == TEST_NUM)
|
|
fprintf(stderr, "\tsupposed answer [%d]: %lf\n", i, yref[i]);
|
|
}
|
|
double ratio = ((double) local_nnz) / ((double) remote_nnz + local_nnz);
|
|
double rows_remote_ratio = ((double) rows_with_nonlocal) / ((double) this->nRows_p);
|
|
if (this->Get_this_pe() == 0) {
|
|
printf("\nRANK 0: global NNZ = %lu\n", remote_nnz + local_nnz);
|
|
printf("RANK 0: global Rows = %d\n", nRows);
|
|
}
|
|
printf("\nLOCALITY %d : Remote/Local cols %lu/%lu Fraction Columns Local %f Fraction Rows with Remote Columns %f\n", this->Get_this_pe(), remote_nnz, local_nnz, ratio, rows_remote_ratio);
|
|
}
|
|
|
|
template<typename FloatType>
|
|
bool SparseTriangularSolve<FloatType>::CSRCheckCPU(FloatType alpha)
|
|
{
|
|
FloatType *NZvalues = SparseMatrix<FloatType>::GetVals();
|
|
int *Cols = SparseMatrix<FloatType>::GetCols();
|
|
int *rowptrs = SparseMatrix<FloatType>::GetRowPtrs();
|
|
double internal_alpha = alpha;
|
|
bool all_worked = true;
|
|
|
|
#pragma omp parallel for
|
|
for(int i = 0; i < nRows; i++)
|
|
{
|
|
#pragma omp flush (all_worked)
|
|
if (all_worked)
|
|
{
|
|
double temp = 0.;
|
|
for(int j = rowptrs[i]; j < rowptrs[i+1]; j++)
|
|
{
|
|
int ci = Cols[j];
|
|
// Skip anything that lies on the diagonal. We need to solve for that.
|
|
temp += NZvalues[j] * yref[ci];
|
|
}
|
|
double compare_val = 0.;
|
|
double alpha_x = internal_alpha * x[i];
|
|
if(typeid(FloatType) == typeid(float))
|
|
{
|
|
compare_val = fabs(alpha_x*1e-3);
|
|
if (compare_val < 10*FLT_EPSILON)
|
|
compare_val = 10*FLT_EPSILON;
|
|
if ((FloatType)(alpha_x - compare_val) > (FloatType)temp || (FloatType)(alpha_x + compare_val) < (FloatType)temp)
|
|
{
|
|
fflush(stdout);
|
|
fprintf(stderr, " CPU CALCULATION ERROR on row %d\n", i);
|
|
fprintf(stderr, "\tReal value for row %d: %.*e\n", i, OP_FLT_Digs-1, (float)alpha_x);
|
|
fprintf(stderr, "\tCalculated value for row %d: %.*e\n", i, OP_FLT_Digs-1, (float)temp);
|
|
all_worked = false;
|
|
#pragma omp flush (all_worked)
|
|
}
|
|
}
|
|
else if(typeid(FloatType) == typeid(double))
|
|
{
|
|
compare_val = fabs(alpha_x*1e-4);
|
|
if (compare_val < 10*DBL_EPSILON)
|
|
compare_val = 10*DBL_EPSILON;
|
|
if ((FloatType)(alpha_x - compare_val) > (FloatType)temp || (FloatType)(alpha_x + compare_val) < (FloatType)temp)
|
|
{
|
|
fflush(stdout);
|
|
fprintf(stderr, " CPU CALCULATION ERROR on row %d\n", i);
|
|
fprintf(stderr, "\tReal value for row %d: %.*le\n", i, OP_DBL_Digs-1, (double)alpha_x);
|
|
fprintf(stderr, "\tCalculated value for row %d: %.*le\n", i, OP_DBL_Digs-1, (double)temp);
|
|
all_worked = false;
|
|
#pragma omp flush (all_worked)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return all_worked;
|
|
}
|
|
|
|
template<>
|
|
int SparseTriangularSolve<float>::VerifyResults(int iteration)
|
|
{
|
|
int errors = 0;
|
|
|
|
#pragma omp parallel for
|
|
for (int i = 0; i < nRows; i++)
|
|
{
|
|
int assigned_pe = (i / SPTS_BLOCK_SIZE) % this->Get_total_pes();
|
|
if (this->Get_this_pe() == assigned_pe) {
|
|
float compare_val = fabs(yref[i]*1e-3);
|
|
if (compare_val < 10*FLT_EPSILON)
|
|
compare_val = 10*FLT_EPSILON;
|
|
if ((yref[i] - compare_val) > y[i] || (yref[i] + compare_val) < y[i])
|
|
{
|
|
#pragma omp critical
|
|
{
|
|
if(errors == 0)
|
|
{
|
|
fflush(stdout);
|
|
fprintf(stderr, "\nDetected some differences between CPU and GPU results on iteration %d...", iteration);
|
|
}
|
|
fprintf(stderr, "%d GPU CALCULATION ERROR on row %d\n", this->Get_this_pe(), i);
|
|
fprintf(stderr, "\tCPU value for y[%d]: %.*e\n", i, OP_FLT_Digs-1, yref[i]);
|
|
fprintf(stderr, "\tGPU value for y[%d]: %.*e\n", i, OP_FLT_Digs-1, y[i]);
|
|
errors += 1;
|
|
observed_errors[iteration].insert(std::pair<int, float> (i, y[i]));
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return errors;
|
|
}
|
|
|
|
template<>
|
|
int SparseTriangularSolve<double>::VerifyResults(int iteration)
|
|
{
|
|
int errors = 0;
|
|
#pragma omp parallel for
|
|
for (int i = 0; i < nRows; i++)
|
|
{
|
|
double compare_val = fabs(yref[i]*1e-4);
|
|
if (compare_val < 10*DBL_EPSILON)
|
|
compare_val = 10*DBL_EPSILON;
|
|
if ((yref[i] - compare_val) > y[i] || (yref[i] + compare_val) < y[i])
|
|
{
|
|
#pragma omp critical
|
|
{
|
|
if(errors == 0)
|
|
{
|
|
fflush(stdout);
|
|
fprintf(stderr, "\nDetected differences between CPU and GPU results on iteration %d...", iteration);
|
|
}
|
|
fprintf(stderr, "GPU CALCULATION ERROR on row %d\n", i);
|
|
fprintf(stderr, "\tCPU value for y[%d]: %.*e\n", i, OP_DBL_Digs-1, yref[i]);
|
|
fprintf(stderr, "\tGPU value for y[%d]: %.*e\n", i, OP_DBL_Digs-1, y[i]);
|
|
errors += 1;
|
|
observed_errors[iteration].insert(std::pair<int, double> (i, y[i]));
|
|
}
|
|
}
|
|
}
|
|
return errors;
|
|
}
|
|
|
|
template<typename FloatType>
|
|
int SparseTriangularSolve<FloatType>::NonDeterministicErrors()
|
|
{
|
|
int iter = InputFlags::GetValueInt("iterations");
|
|
int non_det_errors = 0;
|
|
#ifdef ALL_SYNCFREE
|
|
for (int i = 1; i < iter; i++)
|
|
{
|
|
if (errors_seen[i] != errors_seen[0])
|
|
{
|
|
non_det_errors++;
|
|
if (non_det_errors == 1)
|
|
{
|
|
fprintf(stderr, "Different SpTS iterations saw different error counts -- non-deterministic bug possible.\n");
|
|
fprintf(stderr, "\te.g. saw %d errors during iteration 0. Saw %d errors during iteration %i\n", errors_seen[0], errors_seen[i], i);
|
|
}
|
|
}
|
|
else if (observed_errors[i] != observed_errors[0])
|
|
{
|
|
non_det_errors++;
|
|
if (non_det_errors == 1)
|
|
{
|
|
fprintf(stderr, "ERRORS were seen. Different iterations saw errors on different rows -- non-deterministic bug possible.\n");
|
|
fprintf(stderr, "\te.g. Iterations 0 and %d were different.\n", i);
|
|
}
|
|
}
|
|
}
|
|
#else
|
|
if (iter >= 1)
|
|
{
|
|
if (errors_seen[0] != errors_seen[1])
|
|
{
|
|
non_det_errors++;
|
|
fprintf(stderr, "Different SpTS iterations saw different error counts -- non-deterministic bug possible.\n");
|
|
fprintf(stderr, "\te.g. saw %d errors during iteration 0. Saw %d errors during iteration %i\n", errors_seen[0], errors_seen[1], 1);
|
|
}
|
|
}
|
|
for (int i = 2; i < iter; i++)
|
|
{
|
|
if (errors_seen[i] != errors_seen[1])
|
|
{
|
|
non_det_errors++;
|
|
if (non_det_errors == 1)
|
|
{
|
|
fprintf(stderr, "Different SpTS iterations saw different error counts -- non-deterministic bug possible.\n");
|
|
fprintf(stderr, "\te.g. saw %d errors during iteration 1. Saw %d errors during iteration %i\n", errors_seen[1], errors_seen[i], i);
|
|
}
|
|
}
|
|
else if (observed_errors[i] != observed_errors[1])
|
|
{
|
|
non_det_errors++;
|
|
if (non_det_errors == 1)
|
|
{
|
|
fprintf(stderr, "ERRORS were seen. Different iterations saw errors on different rows -- non-deterministic bug possible.\n");
|
|
fprintf(stderr, "\te.g. Iterations 1 and %d were different.\n", i);
|
|
}
|
|
}
|
|
}
|
|
#endif
|
|
return non_det_errors;
|
|
}
|
|
|
|
template<typename FloatType>
|
|
int SparseTriangularSolve<FloatType>::MaxErrors()
|
|
{
|
|
int iter = InputFlags::GetValueInt("iterations");
|
|
int max_errors = 0;
|
|
for (int i = 0; i < iter; i++)
|
|
{
|
|
if (errors_seen[i] > max_errors)
|
|
max_errors = errors_seen[i];
|
|
}
|
|
return max_errors;
|
|
}
|
|
|
|
#ifndef __has_builtin
|
|
#define __has_builtin(x) 0
|
|
#endif
|
|
|
|
static inline unsigned int flp2(unsigned int x)
|
|
{
|
|
x |= (x >> 1);
|
|
x |= (x >> 2);
|
|
x |= (x >> 4);
|
|
x |= (x >> 8);
|
|
x |= (x >> 16);
|
|
return x - (x >> 1);
|
|
}
|
|
|
|
// Short rows in CSR-Adaptive are batched together into a single row block.
|
|
// If there are a relatively small number of these, then we choose to do
|
|
// a horizontal reduction (groups of threads all reduce the same row).
|
|
// If there are many threads (e.g. more threads than the maximum size
|
|
// of our workgroup) then we choose to have each thread serially reduce
|
|
// the row.
|
|
// This function calculates the number of threads that could team up
|
|
// to reduce these groups of rows. For instance, if you have a
|
|
// workgroup size of 256 and 4 rows, you could have 64 threads
|
|
// working on each row. If you have 5 rows, only 32 threads could
|
|
// reliably work on each row because our reduction assumes power-of-2.
|
|
template< typename rowBlockType >
|
|
static inline rowBlockType numThreadsForReduction(const rowBlockType num_rows)
|
|
{
|
|
#if defined(__INTEL_COMPILER)
|
|
return 256 >> (_bit_scan_reverse(num_rows-1)+1);
|
|
#elif (defined(__clang__) && __has_builtin(__builtin_clz)) || \
|
|
!defined(__clang) && \
|
|
defined(__GNUG__) && ((__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) > 30202)
|
|
return (256 >> (8*sizeof(int)-__builtin_clz(num_rows-1)));
|
|
#elif defined(_MSC_VER) && (_MSC_VER >= 1400)
|
|
unsigned long bit_returned;
|
|
_BitScanReverse(&bit_returned, (num_rows-1));
|
|
return 256 >> (bit_returned+1);
|
|
#else
|
|
return flp2(256/num_rows);
|
|
#endif
|
|
}
|
|
|
|
template<typename FloatType>
|
|
int SparseTriangularSolve<FloatType>::ComputeRowBlocks(std::vector<uint64_t> &rowBlocks,
|
|
int *rowDelimiters,
|
|
int nRows)
|
|
{
|
|
rowBlocks.erase(rowBlocks.begin(), rowBlocks.end());
|
|
rowBlocks.push_back(0);
|
|
uint64_t sum = 0;
|
|
uint64_t i, last_i = 0;
|
|
|
|
// Check to ensure nRows can fit in 32 bits
|
|
if ((uint64_t) nRows > (uint64_t)pow(2, ROW_BITS))
|
|
{
|
|
fflush(stdout);
|
|
fprintf(stderr, "\nNumber of Rows in the Sparse Matrix is greater than what is supported at present (%d bits) !", ROW_BITS );
|
|
exit(0);
|
|
}
|
|
|
|
int consecutive_long_rows = 0;
|
|
for(i = 1; i <= nRows; i++)
|
|
{
|
|
int row_length = ( rowDelimiters[ i ] - rowDelimiters[ i - 1 ] );
|
|
sum += row_length;
|
|
|
|
// The following section of code calculates whether you're moving between
|
|
// a series of "short" rows and a series of "long" rows.
|
|
// This is because the reduction in CSR-Adaptive likes things to be
|
|
// roughly the same length. Long rows can be reduced horizontally.
|
|
// Short rows can be reduced one-thread-per-row. Try not to mix them.
|
|
if ( row_length > 128 )
|
|
consecutive_long_rows++;
|
|
else if ( consecutive_long_rows > 0 )
|
|
{
|
|
// If it turns out we WERE in a long-row region, cut if off now.
|
|
if (row_length < 32) // Now we're in a short-row region
|
|
consecutive_long_rows = -1;
|
|
else
|
|
consecutive_long_rows++;
|
|
}
|
|
|
|
// If you just entered into a "long" row from a series of short rows,
|
|
// then we need to make sure we cut off those short rows. Put them in
|
|
// their own workgroup.
|
|
if ( consecutive_long_rows == 1 )
|
|
{
|
|
// Assuming there *was* a previous workgroup. If not, nothing to do here.
|
|
if( i - last_i > 1 )
|
|
{
|
|
rowBlocks.push_back( (i - 1) << (64 - ROW_BITS) );
|
|
// If this row fits into CSR-Stream, calculate how many rows
|
|
// can be used to do a parallel reduction.
|
|
// Fill in the low-order bits with the numThreadsForRed
|
|
if (((i-1) - last_i) > 2)
|
|
rowBlocks[rowBlocks.size() - 2] |= numThreadsForReduction((i - 1) - last_i);
|
|
|
|
last_i = i-1;
|
|
sum = row_length;
|
|
}
|
|
}
|
|
else if (consecutive_long_rows == -1)
|
|
{
|
|
// We see the first short row after some long ones that
|
|
// didn't previously fill up a row block.
|
|
rowBlocks.push_back( (i - 1) << (64 - ROW_BITS) );
|
|
if (((i-1) - last_i) > 2)
|
|
rowBlocks[rowBlocks.size() - 2] |= numThreadsForReduction((i - 1) - last_i);
|
|
|
|
last_i = i-1;
|
|
sum = row_length;
|
|
consecutive_long_rows = 0;
|
|
}
|
|
|
|
// Now, what's up with this row? What did it do?
|
|
|
|
// exactly one row results in non-zero elements to be greater than blockSize
|
|
// This is csr-vector case; bottom WG_BITS == workgroup ID
|
|
if( ( i - last_i == 1 ) && sum > 1024 )
|
|
{
|
|
int numWGReq = ceil( (double)row_length / (1024) );
|
|
|
|
// Check to ensure #workgroups can fit in WG_BITS bits, if not
|
|
// then the last workgroup will do all the remaining work
|
|
numWGReq = ( numWGReq < (int)pow( 2, WG_BITS ) ) ? numWGReq : (int)pow( 2, WG_BITS );
|
|
|
|
for( int w = 1; w < numWGReq; w++ )
|
|
{
|
|
rowBlocks.push_back((i-1) << ROW_BITS);
|
|
rowBlocks[rowBlocks.size() - 1] |= (uint64_t)w;
|
|
}
|
|
rowBlocks.push_back(i << ROW_BITS);
|
|
|
|
last_i = i;
|
|
sum = 0;
|
|
consecutive_long_rows = 0;
|
|
}
|
|
// more than one row results in non-zero elements to be greater than blockSize
|
|
// This is csr-stream case; bottom WG_BITS = number of parallel reduction threads
|
|
else if( ( i - last_i > 1 ) && sum > 1024 )
|
|
{
|
|
i--; // This row won't fit, so back off one.
|
|
rowBlocks.push_back( i << (64 - ROW_BITS) );
|
|
if ((i - last_i) > 2)
|
|
rowBlocks[rowBlocks.size() - 2] |= numThreadsForReduction(i - last_i);
|
|
last_i = i;
|
|
sum = 0;
|
|
consecutive_long_rows = 0;
|
|
}
|
|
// This is csr-stream case; bottom WG_BITS = number of parallel reduction threads
|
|
else if( sum == 1024 )
|
|
{
|
|
rowBlocks.push_back( i << (64 - ROW_BITS) );
|
|
if ((i - last_i) > 2)
|
|
rowBlocks[rowBlocks.size() - 2] |= numThreadsForReduction(i - last_i);
|
|
last_i = i;
|
|
sum = 0;
|
|
consecutive_long_rows = 0;
|
|
}
|
|
}
|
|
|
|
// If we didn't fill a row block with the last row, make sure we don't lose it.
|
|
if ( (rowBlocks[rowBlocks.size() - 2] >> (64 - ROW_BITS)) != (uint64_t)(nRows) )
|
|
{
|
|
rowBlocks.push_back( (uint64_t)( nRows ) << (64 - ROW_BITS) );
|
|
if ((nRows - last_i) > 2)
|
|
rowBlocks[rowBlocks.size() - 2] |= numThreadsForReduction(i - last_i);
|
|
}
|
|
|
|
return rowBlocks.size();
|
|
}
|
|
|
|
template<typename FloatType>
|
|
float SparseTriangularSolve<FloatType>::CSRSpTSGPU(uint64_t &ns_per_iter, uint64_t &ns_per_analysis_iter, uint64_t &ns_per_syncfree_iter, uint64_t &ns_per_levelset_iter, uint64_t &ns_per_levelsync_iter, FloatType alpha)
|
|
{
|
|
gpuInt status;
|
|
gpuEvent* event_array;
|
|
#ifdef USE_HIP
|
|
hipSetDevice(this->Get_this_pe());
|
|
hipDeviceProp_t props;
|
|
hipGetDeviceProperties(&props, this->Get_this_pe());
|
|
printf("\nPE %d: PCIe BUS ID %d DEV ID %d\n", this->Get_this_pe(), props.pciBusID, props.pciDeviceID);
|
|
event_array = (gpuEvent*)malloc(sizeof(gpuEvent) * 2);
|
|
hipEventCreate(&event_array[0]);
|
|
hipEventCreate(&event_array[1]);
|
|
#else
|
|
event_array = (gpuEvent*)malloc(sizeof(gpuEvent));
|
|
#endif
|
|
size_t global_work_size;
|
|
size_t local_work_size = WF_PER_WG * WF_SIZE;
|
|
|
|
|
|
/*************************** Setup and create buffers ********************/
|
|
/****** Matrix Setup Code ******/
|
|
/* Get the OpenCL buffers for the input matrix */
|
|
memPointer bufNonZeroes = SparseMatrix<FloatType>::GetDevVals();
|
|
memPointer bufColumnIndices = SparseMatrix<FloatType>::GetDevCols();
|
|
memPointer bufRowPtrs = SparseMatrix<FloatType>::GetDevRowPtrs();
|
|
/* Get the host buffers for the input matrix */
|
|
FloatType *Avalues = SparseMatrix<FloatType>::GetVals();
|
|
int *Acols = SparseMatrix<FloatType>::GetCols();
|
|
int *rowptrs = SparseMatrix<FloatType>::GetRowPtrs();
|
|
|
|
|
|
/****** Adaptive RowBlocks Setup Code ******/
|
|
numBlocks = ComputeRowBlocks(rowBlocks, rowptrs, nRows);
|
|
rowBlocksDev = this->GPU->AllocateMem("rowBlocks", numBlocks*sizeof(int64_t), GPU_MEM_READ_WRITE, NULL);
|
|
uint64_t completedRows = 0;
|
|
completedRowsDev = this->GPU->AllocateMem("completedRows", sizeof(uint64_t), GPU_MEM_READ_WRITE|GPU_MEM_USE_HOST_PTR, &completedRows);
|
|
|
|
/****** SpTS Meta-Data Setup Code ******/
|
|
/* Set up the OpenCL buffers for the SpTS meta-data */
|
|
// TODO -- is this +1 in doneArray nRows+1 required? Why?
|
|
#ifdef USE_ROC_SHMEM
|
|
doneArrayDev = roc_shmem_malloc((nRows+1)*sizeof(uint32_t));
|
|
reqUpdateArrayDev = roc_shmem_malloc((nRows+1)*sizeof(uint32_t));
|
|
shadowDoneArrayDev = roc_shmem_malloc((nRows+1)*sizeof(uint32_t));
|
|
#else
|
|
doneArrayDev = this->GPU->AllocateMem("doneArray", (nRows+1)*sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
reqUpdateArrayDev = this->GPU->AllocateMem("reqUpdateArray", (nRows+1)*sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
shadowDoneArrayDev = this->GPU->AllocateMem("shadowDoneArray", (nRows+1)*sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
#endif
|
|
remoteInProgressArrayDev = this->GPU->AllocateMem("remoteInProgressArray", (nRows+1)*sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
numRowsAtLevelDev = this->GPU->AllocateMem("numRowsAtLevel", (nRows)*sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
rowMapDev = this->GPU->AllocateMem("rowMap", (nRows+1)*sizeof(uint32_t), GPU_MEM_READ_ONLY, NULL);
|
|
maxDepthDev = this->GPU->AllocateMem("maxDepth", sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
totalSpinDev = this->GPU->AllocateMem("totalSpin", sizeof(uint64_t), GPU_MEM_READ_WRITE, NULL);
|
|
oneBufDev = this->GPU->AllocateMem("oneBuf", sizeof(uint32_t), GPU_MEM_READ_WRITE, NULL);
|
|
/* Set up the host buffers for the SpTS meta-data */
|
|
uint32_t *doneArray = (uint32_t*)calloc((nRows+1), sizeof(uint32_t));
|
|
uint32_t *numRowsAtLevel = (uint32_t*)calloc(nRows, sizeof(uint32_t));
|
|
uint32_t *rowMap = (uint32_t*)calloc((nRows+1), sizeof(uint32_t));
|
|
uint32_t maxDepth = 0;
|
|
uint64_t totalSpin = 0;
|
|
|
|
uint32_t *nrows_plus1_zero = (uint32_t*)calloc((nRows+1), sizeof(uint32_t));
|
|
uint64_t u64_zero = 0;
|
|
uint32_t u32_zero = 0;
|
|
|
|
//uint32_t uns_int_one = 0x42280000;
|
|
uint32_t u32_one = 1;
|
|
|
|
// TODO: Gather and flatten out Avalues, Acols, and rowptrs based on
|
|
// row cyclic decomposition. For now, we just copy the hole vals, cols,
|
|
// and row_ptrs matrix, even though we really only access 1/num_pes of the
|
|
// whole thing. We can do some more sophisticated stuff here if we run out
|
|
// of space on the GPU or we don't like the copy overheads of the initial
|
|
// buffers.
|
|
|
|
/************************ Copy initial buffers to device *****************/
|
|
/****** Copy matrix ******/
|
|
this->GPU->CopyToDevice(bufNonZeroes, Avalues, this->nNZ*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(bufColumnIndices, Acols, this->nNZ*sizeof(int), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(bufRowPtrs, rowptrs, (this->nRows+1)*sizeof(int), 0, GPU_TRUE, NULL);
|
|
|
|
/****** Copy vectors ******/
|
|
this->GPU->CopyToDevice(xDev, x, nCols*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(yDev, y_zero, nRows*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
|
|
/****** Copy adaptive rowBlock information ******/
|
|
this->GPU->CopyToDevice(rowBlocksDev, rowBlocks.data(), numBlocks*sizeof(int64_t), 0, GPU_TRUE, NULL);
|
|
|
|
/****** Copy SpTS meta-data needed for analyze_and_solve run ******/
|
|
this->GPU->CopyToDevice(doneArrayDev, nrows_plus1_zero, (nRows+1)*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(shadowDoneArrayDev, nrows_plus1_zero, (nRows+1)*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(reqUpdateArrayDev, nrows_plus1_zero, (nRows+1)*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(remoteInProgressArrayDev, nrows_plus1_zero, (nRows+1)*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(numRowsAtLevelDev, nrows_plus1_zero, nRows*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(maxDepthDev, &u32_zero, sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(totalSpinDev, &u64_zero, sizeof(uint64_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->CopyToDevice(oneBufDev, &u32_one, sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
|
|
/************************** Set up iteration printing ********************/
|
|
/* We want to print, ideally, every iteration that gets up 10% closer to
|
|
* completion. This sets that up */
|
|
int iter = InputFlags::GetValueInt("iterations");
|
|
double print_iter = (float)iter / 10.;
|
|
if (print_iter < 1.)
|
|
print_iter = 1.;
|
|
double next_to_print = 0.;
|
|
|
|
|
|
/**************************** Set up perf analysis ************************/
|
|
// For performance analysis, keep track of how much time we've spent doing
|
|
// kernel work.
|
|
// TODO -- Take in from the command line whether to get kernel or total time.
|
|
// If doing total time, try launching all of the kernels at once and waiting
|
|
// outside. This will apparently reduce the overheads.
|
|
uint64_t total_kern_time = 0;
|
|
uint64_t analyze_kern_time = 0;
|
|
double analyze_kern_flops = 0.;
|
|
uint64_t syncfree_kern_time = 0;
|
|
uint64_t levelset_kern_time = 0;
|
|
uint64_t levelsync_kern_time = 0;
|
|
|
|
errors_seen = new int[iter];
|
|
|
|
int analysis_iter = 0;
|
|
int syncfree_iter = 0;
|
|
int levelset_iter = 0;
|
|
int levelsync_iter = 0;
|
|
|
|
int level_sync_cutoff = InputFlags::GetValueInt("levelsync_size");
|
|
bool syncfree_better = false;
|
|
|
|
int total_workitems_per_workgroup = WF_SIZE * WF_PER_WG;
|
|
//bool roc_shmem_initialized = false;
|
|
|
|
/*********************** Actual work of the benchmark *********************/
|
|
for(int i = 0; i < iter; i++)
|
|
{
|
|
if (i == (int)next_to_print || i == (iter - 1))
|
|
{
|
|
printf("%d..", i+1);fflush(stdout);
|
|
next_to_print += print_iter;
|
|
}
|
|
|
|
#ifndef ALL_SYNCFREE
|
|
#ifdef ALL_ANALYZE
|
|
// When we only want to run the analyze-and-solve mechanism, rather than
|
|
// the more optimized syncfree algorithm, we always go into here.
|
|
if (1)
|
|
#else
|
|
// In any version of the program that has the possibility of running the
|
|
// level-set algorithm, we need to start with the syncfree-and-analyze
|
|
// version of the program, so that we can set up the potential to run the
|
|
// level-set algorithm. This will take place on the first iteration.
|
|
if (i == 0)
|
|
#endif
|
|
{
|
|
analysis_iter++;
|
|
global_work_size = nRows * WF_SIZE;
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_analyze, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
doneArrayDev,
|
|
numRowsAtLevelDev,
|
|
maxDepthDev,
|
|
totalSpinDev);
|
|
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_analyze, 1, NULL, &global_work_size, NULL, 0, NULL, &event_array[0]);
|
|
CL->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
this->GPU->Flush();
|
|
total_kern_time += CL->ComputeTime(event_array[0]);
|
|
analyze_kern_time += CL->ComputeTime(event_array[0]);
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
#ifdef USE_ROC_SHMEM
|
|
global_work_size = this->nRows_p * WF_SIZE;
|
|
num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
/*
|
|
int roc_shmem_queues = (2560 / WF_PER_WG);
|
|
if (2560 % WF_PER_WG)
|
|
roc_shmem_queues++;
|
|
if (!roc_shmem_initialized) {
|
|
int num_threads = InputFlags::GetValueInt("num_roshmem_threads");
|
|
roc_shmem_init(&handle, roc_shmem_queues);
|
|
roc_shmem_initialized = true;
|
|
}
|
|
*/
|
|
int rocshmem_algorithm = InputFlags::GetValueInt("rocshmem_algorithm");
|
|
int rocshmem_put_block_size = InputFlags::GetValueInt("put_block_size");
|
|
int rocshmem_get_backoff_factor = InputFlags::GetValueInt("get_backoff_factor");
|
|
switch (rocshmem_algorithm) {
|
|
case 0:
|
|
printf("Using Put-based intra-kernel algorithm\n");
|
|
break;
|
|
case 1:
|
|
printf("Using Get-based intra-kernel algorithm (Backoff factor %d)\n", rocshmem_get_backoff_factor);
|
|
break;
|
|
case 2:
|
|
printf("Using blocked Put-based intra-kernel algorithm\n");
|
|
printf("Using blocked Put-based intra-kernel algorithm (Block Size %d)\n", rocshmem_put_block_size);
|
|
break;
|
|
case 3:
|
|
printf("Using put/get hybrid intra-kernel algorithm\n");
|
|
break;
|
|
default:
|
|
printf("Unknown ROC_SHMEM algoirthm\n");
|
|
exit(-1);
|
|
}
|
|
size_t LDS_size;
|
|
roc_shmem_dynamic_shared(&LDS_size);
|
|
printf("Work size %zu, wg size %d num workgroups %d LDS %zu thisPE %d Global %d \n", global_work_size, total_workitems_per_workgroup, num_of_workgroups, LDS_size, this->Get_this_pe(), this->Get_total_pes());
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
hipEventRecord(event_array[0], NULL);
|
|
hipLaunchKernelGGL(amd_spts_analyze_and_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
LDS_size, 0,
|
|
global_work_size,
|
|
this->Get_this_pe(),
|
|
this->Get_total_pes(),
|
|
static_cast<unsigned int *>(shadowDoneArrayDev),
|
|
static_cast<unsigned int *>(reqUpdateArrayDev),
|
|
static_cast<unsigned int *>(remoteInProgressArrayDev),
|
|
static_cast<unsigned int *>(oneBufDev),
|
|
rocshmem_algorithm,
|
|
rocshmem_put_block_size,
|
|
rocshmem_get_backoff_factor,
|
|
SPTS_BLOCK_SIZE,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(doneArrayDev),
|
|
static_cast<unsigned int *>(numRowsAtLevelDev),
|
|
static_cast<unsigned int *>(maxDepthDev),
|
|
static_cast<unsigned long long *>(totalSpinDev));
|
|
#else
|
|
hipEventRecord(event_array[0], NULL);
|
|
hipLaunchKernelGGL(amd_spts_analyze_and_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(doneArrayDev),
|
|
static_cast<unsigned int *>(numRowsAtLevelDev),
|
|
static_cast<unsigned int *>(maxDepthDev),
|
|
static_cast<unsigned long long *>(totalSpinDev));
|
|
#endif
|
|
hipEventRecord(event_array[1], NULL);
|
|
hipEventSynchronize(event_array[1]);
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
// Wait for any outstanding network messages to finish up. We
|
|
// can have straggler updates to the doneArray that we don't
|
|
// have any dependencies for but we still eed it to finish so
|
|
// the below statistics can work correctly.
|
|
//ro_shmem_dump_stats(handle);
|
|
//ro_shmem_reset_stats(handle);
|
|
//sleep(10);
|
|
/* if( this->Get_this_pe() == 0 && (this->Get_total_pes() > 1)){
|
|
PRINT_SQ(get_rtn_handle(handle), 0, 1, 0);
|
|
PRINT_CQ(get_rtn_handle(handle), 0, 1, 0);
|
|
PRINT_SQ(get_rtn_handle(handle), 0, 1, 1);
|
|
PRINT_CQ(get_rtn_handle(handle), 0, 1, 1);
|
|
PRINT_SQ(get_rtn_handle(handle), 0, 1, 2);
|
|
PRINT_CQ(get_rtn_handle(handle), 0, 1, 2);
|
|
}*/
|
|
MPI_Barrier(MPI_COMM_WORLD);
|
|
#endif
|
|
float elapsed;
|
|
hipEventElapsedTime(&elapsed, event_array[0], event_array[1]);
|
|
total_kern_time += elapsed * 1000000;
|
|
analyze_kern_time += elapsed * 1000000;
|
|
#endif
|
|
analyze_kern_flops = (2 * (double)nNZ * 1000000000.) / (double)analyze_kern_time;
|
|
this->GPU->CopyToHost(yDev, y, nRows*sizeof(FloatType), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToHost(maxDepthDev, &maxDepth, sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToHost(doneArrayDev, doneArray, (nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToHost(totalSpinDev, &totalSpin, sizeof(uint64_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToHost(numRowsAtLevelDev, numRowsAtLevel, nRows*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
this->GPU->Flush();
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
// Combine global statistics
|
|
MPI_Allreduce(MPI_IN_PLACE, (void *) &maxDepth, 1, MPI_UNSIGNED, MPI_MAX, MPI_COMM_WORLD);
|
|
MPI_Allreduce(MPI_IN_PLACE, (void *) &totalSpin, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
|
|
|
|
// TODO: Broadcast out the doneArray and yDev values to all nodes. This
|
|
// is needed for the below calculations to work since in the 'pull'
|
|
// distributed model we don't request data for rows that we don't
|
|
// have a dependency on.
|
|
#endif
|
|
|
|
bool verify = InputFlags::GetValueBool("verify");
|
|
if (verify) {
|
|
printf("Performing results verification\n");
|
|
errors_seen[i] = VerifyResults(i);
|
|
}
|
|
printf("\nTotalSpin: %lu\n", totalSpin);
|
|
|
|
/* Prefix sum of the number of rows at each level, so that we can
|
|
* calculate how much to offset each level into the rowMap */
|
|
// TODO -- Do this prefix sum on the GPU while copying maxDepth and
|
|
// doneArray back into the host. Set non-blocking on the previous ones.
|
|
this->GPU->CopyToHost(numRowsAtLevelDev, numRowsAtLevel, nRows*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
for (unsigned int joe = 1; joe < maxDepth; joe++)
|
|
numRowsAtLevel[joe] = numRowsAtLevel[joe] + numRowsAtLevel[joe-1];
|
|
|
|
/* Build up the rowMap so that each iteration of the no-wait solve
|
|
* knows what it's global_id->row mapping is.
|
|
* The general mechanism for this is as follows:
|
|
* doneArray[row] holds the level that a particular row is in.
|
|
*
|
|
* We know the total number of levels needed (maxDepth), so rowMap
|
|
* has maxDepth 'buckets'.
|
|
*
|
|
* The numRowsAtLevel array (after the above prefix-sum) tells us
|
|
* how many values are in all of the previous buckets, so that
|
|
* we can get an appropriate array offset for each bucket.
|
|
*
|
|
* The counters array keeps track of how many items are in each
|
|
* bucket so far. Add this to the numRowsAtLevel[] offset.
|
|
*
|
|
* As we walk through all the rows, we check to see which level's
|
|
* bucket we should put this row in. Add it at the end of the
|
|
* current bucket, then increment the counter. */
|
|
uint32_t *counters = (uint32_t *)calloc(maxDepth, sizeof(uint32_t));
|
|
/* for (unsigned int this_row = 0; this_row < nRows; this_row++)
|
|
{
|
|
// We must subtract one here, because the first level is '1'
|
|
// The GPU kernel does that because a value of '0' means
|
|
// 'not done, keep waiting' in the analysis kernel.
|
|
assert(doneArray[this_row] != 0);
|
|
unsigned int this_rows_level = doneArray[this_row] - 1;
|
|
unsigned int previous_level = this_rows_level - 1;
|
|
unsigned int depth_offset;
|
|
if (this_rows_level == 0) // can't check previous level
|
|
depth_offset = 0;
|
|
else
|
|
depth_offset = numRowsAtLevel[previous_level];
|
|
rowMap[depth_offset + counters[this_rows_level]] = this_row;
|
|
counters[this_rows_level] += 1;
|
|
} */
|
|
free(counters);
|
|
this->GPU->CopyToDevice(rowMapDev, rowMap, (nRows+1)*sizeof(uint32_t), 0, GPU_TRUE, NULL);
|
|
free(event_array);
|
|
#ifdef USE_HIP
|
|
event_array = (gpuEvent*)malloc(maxDepth * sizeof(gpuEvent) * 2);
|
|
for (int i = 0; i < maxDepth * 2; i++)
|
|
hipEventCreate(&event_array[i]);
|
|
#else
|
|
event_array = (gpuEvent*)malloc(maxDepth * sizeof(gpuEvent));
|
|
#endif
|
|
#ifdef ALL_ANALYZE
|
|
// We will be coming back into this kernel. Time to reset its data.
|
|
if (i != (iter - 1))
|
|
{
|
|
this->GPU->CopyToDevice(maxDepthDev, &u32_zero, sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(totalSpinDev, &u64_zero, sizeof(uint64_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(numRowsAtLevelDev, nrows_plus1_zero, nRows*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
}
|
|
#endif
|
|
this->GPU->CopyToDevice(yDev, y_zero, nRows*sizeof(FloatType), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(doneArrayDev, nrows_plus1_zero,(nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(shadowDoneArrayDev, nrows_plus1_zero,(nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(remoteInProgressArrayDev, nrows_plus1_zero,(nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(reqUpdateArrayDev, nrows_plus1_zero, (nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->Flush();
|
|
// Either we always want to run just this function block
|
|
// (ALL_ANALYZE), or the first iteration is the analyze-and-solve
|
|
// kernel. Either way, don't continue to the code below this time.
|
|
continue;
|
|
}
|
|
#endif
|
|
// If ALL_SYNCFREE is defined, we always run the amd_spts_syncfree_solve
|
|
// kernel. We never try to speed it up by paying attention to the output
|
|
// levels and running the levelset kernel.
|
|
// If ALL_LEVELSET is set, we only run the analysis kernel up above to get the
|
|
// level-set and do the first solve; after that we skip the
|
|
// amd_spts_analyze_and_solve kernel and do the level-set based solve.
|
|
// Otherwise, we dynamically choose between those kernels based on some
|
|
// statistics that we gathered during the analyze-and-solve run.
|
|
#ifdef ALL_SYNCFREE
|
|
if (1) // always run syncfree algorithm
|
|
#elif defined(ALL_LEVELSET) || defined(ALL_LEVELSYNC)
|
|
if (0) // always *do not* run syncfree algorithm
|
|
#else
|
|
if (totalSpin == 0 || analyze_kern_flops/totalSpin > 25000 || syncfree_better) // Try to run syncfree
|
|
#endif
|
|
{
|
|
syncfree_iter++;
|
|
// TODO -- Eventually get this working with numer of RowBlocks - 1
|
|
global_work_size = nRows * WF_SIZE;
|
|
|
|
uint32_t current_iteration = 0;
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
fprintf(stderr, "ROC_SHMEM not supported for selected algorithm\n");
|
|
exit(-1);
|
|
#endif
|
|
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
doneArrayDev,
|
|
numRowsAtLevelDev,
|
|
maxDepthDev,
|
|
totalSpinDev);
|
|
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel, 1, NULL, &global_work_size, NULL, 0, NULL, &event_array[0]);
|
|
CL->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
current_iteration++;
|
|
this->GPU->Flush();
|
|
total_kern_time += CL->ComputeTime(event_array[0]);
|
|
syncfree_kern_time += CL->ComputeTime(event_array[0]);
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[0], NULL);
|
|
hipLaunchKernelGGL(amd_spts_syncfree_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(doneArrayDev),
|
|
static_cast<unsigned int *>(numRowsAtLevelDev),
|
|
static_cast<unsigned int *>(maxDepthDev),
|
|
static_cast<unsigned long long *>(totalSpinDev));
|
|
hipEventRecord(event_array[1], NULL);
|
|
hipEventSynchronize(event_array[1]);
|
|
current_iteration++;
|
|
float elapsed;
|
|
hipEventElapsedTime(&elapsed, event_array[0], event_array[1]);
|
|
total_kern_time += elapsed * 1000000;
|
|
syncfree_kern_time += elapsed * 1000000;
|
|
|
|
#endif
|
|
|
|
this->GPU->CopyToHost(yDev, y, nRows*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
errors_seen[i] = VerifyResults(i);
|
|
|
|
this->GPU->Flush();
|
|
current_iteration = 0;
|
|
completedRows = 0;
|
|
}
|
|
#if defined(ALL_LEVELSYNC)
|
|
else if (1) // always run levelset+syncfree combination
|
|
#elif defined (ALL_LEVELSET)
|
|
else if (0) // Fall through to level-set
|
|
#else
|
|
else if (1) // always run levelset+syncfree, never fall through to level-set only
|
|
#endif
|
|
{
|
|
// This is the "level-sync" algorithm, where we take the level-set
|
|
// information and launch kernels that combine multiple levels
|
|
// together. This allows us to find parallelism to run on the GPU,
|
|
// even if technically there are some data dependencies between the
|
|
// levels. Within the kernel, we use the synchronization-free algorithm
|
|
// to ensure that we get the right answer.
|
|
// This algorithm reduces the spin-loop overhead of the sync-free
|
|
// algorithm if there are many levels, but it finds more parallelism
|
|
// than the pure level-set algorithms which can only run on one CU.
|
|
levelsync_iter++;
|
|
|
|
// Keep track of total kernels we launch so we can watch for events.
|
|
int total_enqueues = 0;
|
|
|
|
/* The rowMap tells each workgroup within the kernel what
|
|
* rows it is working on. However, each each kernel invocation
|
|
* is working on a different level. Each level is in a separate
|
|
* 'bucket' in the rowMap. We must tell each invocation how far
|
|
* into the rowMap it much index. That's the depth_offset.
|
|
* numRowsAtLevel (after the above prefix-sum) tells us how
|
|
* many rows were in all previous levels combined. */
|
|
unsigned int depth_offset = 0;
|
|
unsigned int running_total = 0; // How many rows in this launch
|
|
|
|
if (level_sync_cutoff == 0)
|
|
{
|
|
if (nRows/maxDepth < 32)
|
|
level_sync_cutoff = 2560;
|
|
else
|
|
level_sync_cutoff = 81920;
|
|
}
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
fprintf(stderr, "ROC_SHMEM not supported for selected algorithm\n");
|
|
exit(-1);
|
|
#endif
|
|
|
|
for (int this_level = 0; this_level < maxDepth; this_level++)
|
|
{
|
|
if (this_level != 0 && running_total == 0)
|
|
depth_offset = numRowsAtLevel[this_level-1];
|
|
|
|
running_total = numRowsAtLevel[this_level] - depth_offset;
|
|
|
|
if (running_total >= level_sync_cutoff)
|
|
{
|
|
global_work_size = (running_total + (running_total % WF_PER_WG)) * WF_SIZE;
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_levelsync, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
doneArrayDev,
|
|
rowMapDev,
|
|
depth_offset);
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_levelsync, 1, NULL, &global_work_size, &local_work_size, 0, NULL, &event_array[total_enqueues]);
|
|
this->GPU->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[total_enqueues * 2], NULL);
|
|
hipLaunchKernelGGL(amd_spts_levelsync_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(doneArrayDev),
|
|
static_cast<unsigned int*>(rowMapDev),
|
|
depth_offset);
|
|
hipEventRecord(event_array[total_enqueues * 2 + 1], NULL);
|
|
#endif
|
|
total_enqueues++;
|
|
running_total = 0;
|
|
}
|
|
}
|
|
if (running_total)
|
|
{
|
|
global_work_size = (running_total + (running_total % WF_PER_WG)) * WF_SIZE;
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_levelsync, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
doneArrayDev,
|
|
rowMapDev,
|
|
depth_offset);
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_levelsync, 1, NULL, &global_work_size, &local_work_size, 0, NULL, &event_array[total_enqueues]);
|
|
this->GPU->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[total_enqueues * 2], NULL);
|
|
hipLaunchKernelGGL(amd_spts_levelsync_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(doneArrayDev),
|
|
static_cast<unsigned int*>(rowMapDev),
|
|
depth_offset);
|
|
hipEventRecord(event_array[total_enqueues * 2 + 1], NULL);
|
|
#endif
|
|
total_enqueues++;
|
|
}
|
|
|
|
// After we cross this clFinish, all of the kernel invocations have
|
|
// completed, and the final answer is in yDev. Now we should add up
|
|
// all of the kernel runtimes from all levels to see how long this
|
|
// levelset solve took.
|
|
this->GPU->Flush();
|
|
for (int this_enqueue = 0; this_enqueue < total_enqueues; this_enqueue++)
|
|
{
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
total_kern_time += CL->ComputeTime(event_array[this_enqueue]);
|
|
levelsync_kern_time += CL->ComputeTime(event_array[this_enqueue]);
|
|
#else
|
|
float elapsed;
|
|
hipEventElapsedTime(&elapsed, event_array[this_enqueue * 2], event_array[this_enqueue * 2 + 1]);
|
|
total_kern_time += elapsed * 1000000;
|
|
levelsync_kern_time += elapsed * 1000000;
|
|
#endif
|
|
}
|
|
// The analyze kernel is about 15% slower than the syncfree kernel.
|
|
// As such, if the level-sync verseion is < 15% faster, it's likely
|
|
// that syncfree will win. Let's go back to doing that.
|
|
if (i == 1 && (analyze_kern_time < (levelsync_kern_time * 1.15)))
|
|
syncfree_better = true;
|
|
this->GPU->CopyToHost(yDev, y, nRows*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
errors_seen[i] = VerifyResults(i);
|
|
}
|
|
else // Run level-set algorithm
|
|
{
|
|
// This is the level-set SpTS kernel, which can be done after the
|
|
// first analyze-and-solve kernel. In this case, we know the levels
|
|
// that each row is in, so we can launch one kernel per level with
|
|
// exactly the right number of workgroups (one WG per row).
|
|
// This means that we don't have any in-kernel atomics, spin-loops,
|
|
// etc, so things run much faster. However, we much launch a
|
|
// potentially large number of kernels.
|
|
// Number of levels is maxDepth. */
|
|
levelset_iter++;
|
|
|
|
#ifdef USE_ROC_SHMEM
|
|
fprintf(stderr, "ROC_SHMEM not supported for selected algorithm\n");
|
|
exit(-1);
|
|
#endif
|
|
|
|
// Keep track of total kernels we launch so we can watch for events.
|
|
int total_enqueues = 0;
|
|
|
|
unsigned int start_level = 0;
|
|
unsigned int end_level = 0;
|
|
unsigned int in_a_run = 0;
|
|
unsigned int running_total = 0;
|
|
|
|
// How far into the rowMap that lists which rows are in each level
|
|
unsigned int depth_offset = 0;
|
|
|
|
unsigned int total_vector = 0;
|
|
unsigned int total_levelset = 0;
|
|
for (int this_level = 0; this_level < maxDepth; this_level++)
|
|
{
|
|
unsigned int inner_depth_offset;
|
|
if (this_level == 0)
|
|
inner_depth_offset = 0;
|
|
else
|
|
inner_depth_offset = numRowsAtLevel[this_level-1];
|
|
unsigned int total_in_this_depth = numRowsAtLevel[this_level] - inner_depth_offset;
|
|
|
|
if (total_in_this_depth == 0)
|
|
continue;
|
|
|
|
end_level = this_level;
|
|
// Comment out this if(){} section to force us to always
|
|
// launch the levelset kernel.
|
|
if (total_in_this_depth <= 2*WF_PER_WG)
|
|
{
|
|
running_total += total_in_this_depth;
|
|
if (in_a_run == 0)
|
|
{
|
|
start_level = this_level;
|
|
depth_offset = inner_depth_offset;
|
|
in_a_run = 1;
|
|
}
|
|
}
|
|
else
|
|
{
|
|
if (in_a_run)
|
|
{
|
|
global_work_size = WF_SIZE * WF_PER_WG;
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_vector, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
rowMapDev,
|
|
numRowsAtLevelDev,
|
|
depth_offset,
|
|
start_level,
|
|
end_level);
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_vector, 1, NULL, &global_work_size, &global_work_size, 0, NULL, &event_array[total_enqueues]);
|
|
this->GPU->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[total_enqueues * 2], NULL);
|
|
hipLaunchKernelGGL(amd_spts_vector_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(rowMapDev),
|
|
static_cast<unsigned int *>(numRowsAtLevelDev),
|
|
depth_offset,
|
|
start_level,
|
|
end_level);
|
|
hipEventRecord(event_array[total_enqueues * 2 + 1], NULL);
|
|
#endif
|
|
total_enqueues++;
|
|
//printf("\n\tVector. offset %u Start %u End %u Rows in this enq %u\n", depth_offset, start_level, end_level, running_total);
|
|
in_a_run = start_level = end_level = running_total = 0;
|
|
depth_offset = numRowsAtLevel[this_level-1];
|
|
total_vector++;
|
|
}
|
|
global_work_size = WF_SIZE * total_in_this_depth;
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_levelset, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
rowMapDev,
|
|
depth_offset,
|
|
alpha);
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_levelset, 1, NULL, &global_work_size, NULL, 0, NULL, &event_array[total_enqueues]);
|
|
this->GPU->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[total_enqueues * 2], NULL);
|
|
hipLaunchKernelGGL(amd_spts_levelset_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
static_cast<unsigned int *>(rowMapDev),
|
|
depth_offset,
|
|
alpha);
|
|
hipEventRecord(event_array[total_enqueues * 2 + 1], NULL);
|
|
#endif
|
|
total_enqueues++;
|
|
depth_offset = numRowsAtLevel[this_level];
|
|
total_levelset++;
|
|
}
|
|
}
|
|
end_level++;
|
|
if (in_a_run)
|
|
{
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
CL->SetArgs(CLHelper::SpTSKernel_vector, 0,
|
|
bufNonZeroes,
|
|
bufColumnIndices,
|
|
bufRowPtrs,
|
|
xDev,
|
|
yDev,
|
|
alpha,
|
|
rowMapDev,
|
|
numRowsAtLevelDev,
|
|
depth_offset,
|
|
start_level,
|
|
end_level);
|
|
global_work_size = WF_SIZE * WF_PER_WG;
|
|
status = clEnqueueNDRangeKernel(CLHelper::commandQueue, CLHelper::SpTSKernel_vector, 1, NULL, &global_work_size, &global_work_size, 0, NULL, &event_array[total_enqueues]);
|
|
this->GPU->checkStatus(status,"clEnqueueNDRangeKernel failed");
|
|
#else
|
|
int num_of_workgroups = (global_work_size + total_workitems_per_workgroup - 1)
|
|
/ total_workitems_per_workgroup;
|
|
hipEventRecord(event_array[total_enqueues * 2], NULL);
|
|
hipLaunchKernelGGL(amd_spts_vector_solve,
|
|
dim3(num_of_workgroups),
|
|
dim3(total_workitems_per_workgroup),
|
|
0, 0,
|
|
global_work_size,
|
|
static_cast<FPTYPE *>(bufNonZeroes),
|
|
static_cast<int *>(bufColumnIndices),
|
|
static_cast<int *>(bufRowPtrs),
|
|
static_cast<FPTYPE *>(xDev),
|
|
static_cast<FPTYPE *>(yDev),
|
|
alpha,
|
|
static_cast<unsigned int *>(rowMapDev),
|
|
static_cast<unsigned int *>(numRowsAtLevelDev),
|
|
depth_offset,
|
|
start_level,
|
|
end_level);
|
|
hipEventRecord(event_array[total_enqueues * 2 + 1], NULL);
|
|
#endif
|
|
total_enqueues++;
|
|
//printf("\n\tVector. offset %u Start %u End %u Rows in this enq %u\n", depth_offset, start_level, end_level, running_total);
|
|
in_a_run = start_level = end_level = running_total = 0;
|
|
total_vector++;
|
|
}
|
|
|
|
if (i == 1)
|
|
printf("\nTotal Vector: %u\nTotal levelset: %u\n", total_vector, total_levelset);
|
|
// After we cross this clFinish, all of the kernel invocations have
|
|
// completed, and the final answer is in yDev. Now we should add up
|
|
// all of the kernel runtimes from all levels to see how long this
|
|
// levelset solve took.
|
|
this->GPU->Flush();
|
|
for (int this_enqueue = 0; this_enqueue < total_enqueues; this_enqueue++)
|
|
{
|
|
#ifndef USE_HIP
|
|
CLHelper *CL = dynamic_cast<CLHelper*>(this->GPU);
|
|
total_kern_time += CL->ComputeTime(event_array[this_enqueue]);
|
|
levelset_kern_time += CL->ComputeTime(event_array[this_enqueue]);
|
|
#else
|
|
float elapsed;
|
|
hipEventElapsedTime(&elapsed, event_array[this_enqueue * 2], event_array[this_enqueue * 2 + 1]);
|
|
total_kern_time += elapsed * 1000000;
|
|
levelset_kern_time += elapsed * 1000000;
|
|
#endif
|
|
}
|
|
this->GPU->CopyToHost(yDev, y, nRows*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
errors_seen[i] = VerifyResults(i);
|
|
}
|
|
|
|
#ifndef ALL_SYNCFREE
|
|
if (i == 1)
|
|
printf("\nmaxDepth %d\n", maxDepth);
|
|
#endif
|
|
if (i != (iter - 1))
|
|
{
|
|
this->GPU->CopyToDevice(yDev, y_zero, nRows*sizeof(FloatType), 0, GPU_FALSE, NULL);
|
|
this->GPU->CopyToDevice(doneArrayDev, nrows_plus1_zero,(nRows+1)*sizeof(uint32_t), 0, GPU_FALSE, NULL);
|
|
this->GPU->Flush();
|
|
}
|
|
}
|
|
|
|
float gflops = 0.f;
|
|
printf("\n\nnnz: %d\n", nNZ);
|
|
gflops = (float)(2 * nNZ) / (float)(total_kern_time/iter);
|
|
ns_per_iter = total_kern_time/iter;
|
|
|
|
if (analysis_iter > 0)
|
|
ns_per_analysis_iter = analyze_kern_time / analysis_iter;
|
|
else
|
|
ns_per_analysis_iter = 0;
|
|
if (syncfree_iter > 0)
|
|
ns_per_syncfree_iter = syncfree_kern_time / syncfree_iter;
|
|
else
|
|
ns_per_syncfree_iter = 0;
|
|
if (levelset_iter > 0)
|
|
ns_per_levelset_iter = levelset_kern_time / levelset_iter;
|
|
else
|
|
ns_per_levelset_iter = 0;
|
|
if (levelsync_iter > 0)
|
|
ns_per_levelsync_iter = levelsync_kern_time / levelsync_iter;
|
|
else
|
|
ns_per_levelsync_iter = 0;
|
|
|
|
this->GPU->CopyToHost(yDev, y, nRows*sizeof(FloatType), 0, GPU_TRUE, NULL);
|
|
|
|
if (doneArray)
|
|
free(doneArray);
|
|
if (numRowsAtLevel)
|
|
free(numRowsAtLevel);
|
|
if (rowMap)
|
|
free(rowMap);
|
|
if (event_array)
|
|
free(event_array);
|
|
|
|
return gflops;
|
|
}
|
|
|
|
#endif //SpTS_H
|