- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在构建一个 CUDA 内核来计算数值 N*N
函数的雅可比,使用有限差分;在我提供的示例中,它是平方函数(向量的每个条目都平方)。主机编码在线性内存中分配,而我在内核中使用二维索引。
我的问题是 我还没有找到对角线求和的方法 矩阵 cudaMalloc
'ed。我一直在尝试使用声明 threadIdx.x == blockIdx.x
作为对角线的条件,但它的计算结果为 true
只为他们俩在 0
.
这是内核和 编辑:我根据评论中的建议发布了整个代码作为答案(main()
基本相同,而内核不是)
template <typename T>
__global__ void jacobian_kernel (
T * J,
const T t0,
const T tn,
const T h,
const T * u0,
const T * un,
const T * un_old)
{
T cgamma = 2 - sqrtf(2);
const unsigned int t = threadIdx.x;
const unsigned int b = blockIdx.x;
const unsigned int tid = t + b * blockDim.x;
/*__shared__*/ T temp_sx[BLOCK_SIZE][BLOCK_SIZE];
/*__shared__*/ T temp_dx[BLOCK_SIZE][BLOCK_SIZE];
__shared__ T sm_temp_du[BLOCK_SIZE];
T* temp_du = &sm_temp_du[0];
if (tid < N )
{
temp_sx[b][t] = un[t];
temp_dx[b][t] = un[t];
if ( t == b )
{
if ( tn == t0 )
{
temp_du[t] = u0[t]*0.001;
temp_sx[b][t] += temp_du[t]; //(*)
temp_dx[b][t] -= temp_du[t];
temp_sx[b][t] += ( abs( temp_sx[b][t] ) < 10e-6 ? 0.1 : 0 );
temp_dx[b][t] += ( abs( temp_dx[b][t] ) < 10e-6 ? 0.1 : 0 );
temp_sx[b][t] = ( temp_sx[b][t] == 0 ? 0.1 : temp_sx[b][t] );
temp_dx[b][t] = ( temp_dx[b][t] == 0 ? 0.1 : temp_dx[b][t] );
}
else
{
temp_du[t] = MAX( un[t] - un_old[t], 10e-6 );
temp_sx[b][t] += temp_du[t];
temp_dx[b][t] -= temp_du[t];
}
}
__syncthreads();
//J = f(tn, un + du)
d_func(tn, (temp_sx[b]), (temp_sx[b]), 1.f);
d_func(tn, (temp_dx[b]), (temp_dx[b]), 1.f);
__syncthreads();
J[tid] = (temp_sx[b][t] - temp_dx[b][t]) * powf((2 * temp_du[t]), -1);
//J[tid]*= - h*cgamma/2;
//J[tid]+= ( t == b ? 1 : 0);
//J[tid] = temp_J[tid];
}
}
un
进入temp_sx
的每一行和 temp_dx
du
作为 0.01
震级来自 u0
du
到 temp_sx
的对角线, 减去 du
来自 temp_dx
的对角线temp_sx
的每个条目的平方函数和 temp_dx
2*du
(f(un + du*e_i) - f(un - du*e_i))/2*du
.
du
到矩阵的对角线 的
temp_sx
和
temp_dx
就像我在
(*)
中尝试过的一样.我怎样才能做到这一点?
.y
内核中根本没有使用轴。我正在使用固定数量的共享内存调用内核
int main()
我正在调用内核
#define REAL sizeof(float)
#define N 32
#define BLOCK_SIZE 16
#define NUM_BLOCKS ((N*N + BLOCK_SIZE - 1)/ BLOCK_SIZE)
...
dim3 dimGrid(NUM_BLOCKS,);
dim3 dimBlock(BLOCK_SIZE);
size_t shm_size = N*N*REAL;
jacobian_kernel <<< dimGrid, dimBlock, size_t shm_size >>> (...);
if(threadIdx.x == blockIdx.x){...}
.
为什么这不正确? 我问它是因为在调试和使代码打印语句时,它只计算
true
如果它们都是 0。因此
du[0]
是唯一的数值,矩阵变为
nan
.请注意,这种方法适用于我构建的第一个代码,我在其中调用内核
jacobian_kernel <<< N, N >>> (...)
threadIdx.x == blockIdx.x
元素在对角线上。但是这种方法不再适用,因为现在我需要处理更大的
N
(可能大于 1024,这是每个块的最大线程数)。
最佳答案
根据对答案的评论中的建议,这是我设法解决问题的方法。该示例是可编译的,前提是您输入 helper_cuda.h
和 helper_string.h
在同一目录中或添加 -I
CUDA 示例的指令包含路径,与 CUDA 工具包一起安装。相关的变化只在内核中; main()
有一个小的变化虽然,因为我调用了双倍的资源来执行内核,但是 .y
线程块网格的轴甚至根本没有使用,所以它没有产生任何错误。
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <assert.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include "helper_cuda.h"
#include "helper_string.h"
#include <fstream>
#ifndef MAX
#define MAX(a,b) ((a > b) ? a : b)
#endif
#define REAL sizeof(float)
#define N 128
#define BLOCK_SIZE 128
#define NUM_BLOCKS ((N*N + BLOCK_SIZE - 1)/ BLOCK_SIZE)
template <typename T>
inline void printmatrix( T mat, int rows, int cols);
template <typename T>
__global__ void jacobian_kernel ( const T * A, T * J, const T t0, const T tn, const T h, const T * u0, const T * un, const T * un_old);
template<typename T>
__device__ void d_func(const T t, const T u[], T res[], const T h = 1);
template<typename T>
int main ()
{
float t0 = 0.; //float tn = 0.;
float h = 0.1;
float* u0 = (float*)malloc(REAL*N); for(int i = 0; i < N; ++i){u0[i] = i+1;}
float* un = (float*)malloc(REAL*N); memcpy(un, u0, REAL*N);
float* un_old = (float*)malloc(REAL*N); memcpy(un_old, u0, REAL*N);
float* J = (float*)malloc(REAL*N*N);
float* A = (float*)malloc(REAL*N*N); host_heat_matrix(A);
float *d_u0;
float *d_un;
float *d_un_old;
float *d_J;
float *d_A;
checkCudaErrors(cudaMalloc((void**)&d_u0, REAL*N)); //printf("1: %p\n", d_u0);
checkCudaErrors(cudaMalloc((void**)&d_un, REAL*N)); //printf("2: %p\n", d_un);
checkCudaErrors(cudaMalloc((void**)&d_un_old, REAL*N)); //printf("3: %p\n", d_un_old);
checkCudaErrors(cudaMalloc((void**)&d_J, REAL*N*N)); //printf("4: %p\n", d_J);
checkCudaErrors(cudaMalloc((void**)&d_A, REAL*N*N)); //printf("4: %p\n", d_J);
checkCudaErrors(cudaMemcpy(d_u0, u0, REAL*N, cudaMemcpyHostToDevice)); assert(d_u0 != NULL);
checkCudaErrors(cudaMemcpy(d_un, un, REAL*N, cudaMemcpyHostToDevice)); assert(d_un != NULL);
checkCudaErrors(cudaMemcpy(d_un_old, un_old, REAL*N, cudaMemcpyHostToDevice)); assert(d_un_old != NULL);
checkCudaErrors(cudaMemcpy(d_J, J, REAL*N*N, cudaMemcpyHostToDevice)); assert(d_J != NULL);
checkCudaErrors(cudaMemcpy(d_A, A, REAL*N*N, cudaMemcpyHostToDevice)); assert(d_A != NULL);
dim3 dimGrid(NUM_BLOCKS); std::cout << "NUM_BLOCKS \t" << dimGrid.x << "\n";
dim3 dimBlock(BLOCK_SIZE); std::cout << "BLOCK_SIZE \t" << dimBlock.x << "\n";
size_t shm_size = N*REAL; //std::cout << shm_size << "\n";
//HERE IS A RELEVANT CHANGE OF THE MAIN, SINCE I WAS CALLING
//THE KERNEL WITH A 2D GRID BUT WITHOUT USING THE .y AXIS,
//WHILE NOW THE GRID IS 1D
jacobian_kernel <<< dimGrid, dimBlock, shm_size >>> (d_A, d_J, t0, t0, h, d_u0, d_un, d_un_old);
checkCudaErrors(cudaMemcpy(J, d_J, REAL*N*N, cudaMemcpyDeviceToHost)); //printf("4: %p\n", d_J);
printmatrix( J, N, N);
checkCudaErrors(cudaDeviceReset());
free(u0);
free(un);
free(un_old);
free(J);
}
template <typename T>
__global__ void jacobian_kernel (
const T * A,
T * J,
const T t0,
const T tn,
const T h,
const T * u0,
const T * un,
const T * un_old)
{
T cgamma = 2 - sqrtf(2);
const unsigned int t = threadIdx.x;
const unsigned int b = blockIdx.x;
const unsigned int tid = t + b * blockDim.x;
/*__shared__*/ T temp_sx[BLOCK_SIZE][BLOCK_SIZE];
/*__shared__*/ T temp_dx[BLOCK_SIZE][BLOCK_SIZE];
__shared__ T sm_temp_du;
T* temp_du = &sm_temp_du;
//HERE IS A RELEVANT CHANGE (*)
if ( t < BLOCK_SIZE && b < NUM_BLOCKS )
{
temp_sx[b][t] = un[t]; //printf("temp_sx[%d] = %f\n", t,(temp_sx[b][t]));
temp_dx[b][t] = un[t];
//printf("t = %d, b = %d, t + b * blockDim.x = %d \n",t, b, tid);
//HERE IS A NOTE (**)
if ( t == b )
{
//printf("t = %d, b = %d \n",t, b);
if ( tn == t0 )
{
*temp_du = u0[t]*0.001;
temp_sx[b][t] += *temp_du;
temp_dx[b][t] -= *temp_du;
temp_sx[b][t] += ( abs( temp_sx[b][t] ) < 10e-6 ? 0.1 : 0 );
temp_dx[b][t] += ( abs( temp_dx[b][t] ) < 10e-6 ? 0.1 : 0 );
temp_sx[b][t] = ( temp_sx[b][t] == 0 ? 0.1 : temp_sx[b][t] );
temp_dx[b][t] = ( temp_dx[b][t] == 0 ? 0.1 : temp_dx[b][t] );
}
else
{
*temp_du = MAX( un[t] - un_old[t], 10e-6 );
temp_sx[b][t] += *temp_du;
temp_dx[b][t] -= *temp_du;
}
;
}
//printf("du[%d] %f\n", tid, (*temp_du));
__syncthreads();
//printf("temp_sx[%d][%d] = %f\n", b, t, temp_sx[b][t]);
//printf("temp_dx[%d][%d] = %f\n", b, t, temp_dx[b][t]);
//d_func(tn, (temp_sx[b]), (temp_sx[b]), 1.f);
//d_func(tn, (temp_dx[b]), (temp_dx[b]), 1.f);
matvec_dev( tn, A, (temp_sx[b]), (temp_sx[b]), N, N, 1.f );
matvec_dev( tn, A, (temp_dx[b]), (temp_dx[b]), N, N, 1.f );
__syncthreads();
//printf("temp_sx_later[%d][%d] = %f\n", b, t, (temp_sx[b][t]));
//printf("temp_sx_later[%d][%d] - temp_dx_later[%d][%d] = %f\n", b,t,b,t, (temp_sx[b][t] - temp_dx[b][t]) / 2 * *temp_du);
//if (t == b ) printf( "2du[%d]^-1 = %f\n",t, powf((2 * *temp_du), -1));
J[tid] = (temp_sx[b][t] - temp_dx[b][t]) / (2 * *temp_du);
}
}
template<typename T>
__device__ void d_func(const T t, const T u[], T res[], const T h )
{
__shared__ float temp_u;
temp_u = u[threadIdx.x];
res[threadIdx.x] = h*powf( (temp_u), 2);
}
template <typename T>
inline void printmatrix( T mat, int rows, int cols)
{
std::ofstream matrix_out;
matrix_out.open( "heat_matrix.txt", std::ofstream::out);
for( int i = 0; i < rows; i++)
{
for( int j = 0; j <cols; j++)
{
double next = mat[i + N*j];
matrix_out << ( (next >= 0) ? " " : "") << next << " ";
}
matrix_out << "\n";
}
}
(*)
.在我使用之前
if (tid < N)
这有两个缺点:
tid < N*N
,因为我的数据是二维的,而 tid
是一个跟踪所有数据的全局索引。 tid < N*N
,因为我将函数调用分成块,t < BLOCK_SIZE && b < NUM_BLOCKS
我似乎更清楚索引在代码中的排列方式。 t == b
在
(**)
实际上是对对角元素进行操作的正确方法 的矩阵。它被评估的事实
true
仅在
0
是因为我上面的错误。
关于matrix - 处理 CUDA 中的矩阵 : understanding basic concepts,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/41764331/
以下是 BASIC 中的示例程序。如果标记的条件不为真,有人能告诉我这个函数返回什么吗?我必须将程序移植到 C++ 并且需要理解它。我没有基础知识 - 请忍受简单的问题。 FUNCTION Check
所以,我正在为一个学校项目而苦苦挣扎,我有一个在 PALM 上编程的 BASIC 代码,我需要将其转换为 LabView,所以首先我要将代码转换为伪代码,但我已经以前从未使用过 BASIC,所以我在处
嗯,我突然非常想念 GW-Basic,所以我想在我的新 Ubuntu 盒子里安装一个。我已经很长时间没有为 GW-Basic 编程了,所以我不知道是否有新版本发布。我从旧的 DOS 3.3 盒子中复制
我正在寻找可用于查找 BASIC 语法的引用。该代码用在旧的 HP-85 上,可能会影响命令列表。我在使用谷歌时遇到了问题,因为每次我搜索“<>基本命令”之类的东西时,它都会认为我的字面意思是基本,并
我继承了一个 BASIC 脚本,我正在尝试将其重写为 Python。我不知道 BASIC,甚至不知道这是哪个版本的 BASIC。是 Visual Basic 吗?请帮我翻译这段代码。 'County
我正在编写一个程序,可以判断给定的数字是否是素数。无论我输入素数还是其他数字,总是显示“这不是素数”。这其中有什么错误吗? 10 input "what is the number";a
我继承了一个 BASIC 脚本,我正在尝试将其重写为 Python。我不知道 BASIC,甚至不知道这是哪个版本的 BASIC。是 Visual Basic 吗?请帮我翻译这段代码。 'County
我最近翻出了我的一本旧书,夏威夷计算机之谜,出版于 1985 年。第 81 页有一段 BASIC 代码, 1 For N = 7 to 77 2 Print N, SQR(N) - INT (SQR
在大多数编程语言中,您可以在输出过程中将字符串与变量混合和匹配。但是,我似乎找不到找到这种方法的好方法。这是我的代码: Prompt A,B √(A^2+B^2)->C If iPart(C)≠C T
就目前情况而言,这个问题不太适合我们的问答形式。我们希望答案得到事实、引用资料或专业知识的支持,但这个问题可能会引发辩论、争论、民意调查或扩展讨论。如果您觉得这个问题可以改进并可能重新开放,visit
我正在尝试将用 Rocky Mountain BASIC 编写的程序移植到 GWT,但我被以下语句难住了: 1040 Cfs = 0.75/((LGT(Reyns)-2)^2) 1040是行号。 Re
以下两个(功能相同的)程序摘自旧一期的 Compute's Gazette。主要区别在于程序 1 将目标基本内存位置(7680 和 38400)内联,而程序 2 首先将它们分配给一个变量。 程序 1
每一天都有自己的数据,我需要在不使用最近数据的情况下运行一些科学的东西。 基本上我需要阻止注释掉它。这在 Liberty Basic 中是如何完成的? 最佳答案 Liberty Basic 不支持多行
' Gambas class file ' Math Drill by William Teder. Feel free to use parts of the code, but please gi
我一直在寻找在 TIBASIC 中制作一个程序,该程序可以评估代码正在运行的计算器类型,无需汇编。因为我认为没有任何东西可以从 about 屏幕获取信息。这是我想出的一段代码: :ClrDraw :T
我正在用 TI-Basic 编写一个简单的 Pong 游戏,但编辑器不允许我在我已经编写的代码中插入一行。 例如 print "Hello world" <--Where I want to inse
我是 1980 年代早期/中期个人电脑的忠实粉丝,例如 Amstrad CPC、Commodore 64 和 Sinclair Spectrum。这些计算机都拥有的一件事是 BASIC 版本。 作为一
所以从 70 年代开始就使用 Pick 系统。我们所做的一切都是在 Pick 中完成的。我想维护 Pick 记录,但使用另一种语言(例如 Java)作为前端用户界面。问题是 D3 似乎被锁定在 lin
BASIC 编程语言中 GOTO 和 GOSUB 语句有什么区别? 最佳答案 GOTO 只是跳转到另一行,GOSUB 会跟踪它的来源(大概是在堆栈上),因此当解释器遇到 RETURN 时,它返回到最后
我父亲在 80 年代学习了编程,但他仍然坚持使用 GW-BASIC(并以此谋生)。要求他创建一个 CSV 文件,但他只知道如何创建固定宽度记录的文件。 我在网上发现打开纯文本文件的语法是: OPEN
我是一名优秀的程序员,十分优秀!