c - 如何在 C 和 CUDA 中通过修改后的 Gram-Schmidt 方法进行 QR 分解-6ren

c - 如何在 C 和 CUDA 中通过修改后的 Gram-Schmidt 方法进行 QR 分解

转载作者：太空宇宙更新时间：2023-11-04 01:08:14

<分区>

任何人都知道如何在 C 和 CUDA 中通过改进的 Gram-Schmidt 方法进行 QR 分解。一些例子/来源/论文或其他？非常感谢。

编辑:我无法回答我的问题，因为有人已经关闭了它，所以我决定更新我的问题。

/*
 * QR decomposition via modified Gram-Schmidt algorithm
 * 
 * @Package = QR-decomposition
 * @Program = QR_gpu
 * @Version = 13.0928
 */

// Libraries
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <ctype.h>
#include <time.h>
#include <math.h>
#include <cuda.h>

// Constant
#define PROG "QR_cpu"
#define VERSION "13.1003"
#define PACKAGE "QR-Decomposition"

// Threads per block
#define THREAD_P_BLOCK 512
// Blocks per grid
#define BLOCKS_P_GRID  512

// macro
/* wrap each API call with the gpuErrchk macro, which will process the return 
status of the API call it wraps: http://bit.ly/1dTD0ZE */
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }

// Prototypes
__global__ void xTA (float *, float *, int, int, int, int, int);
__global__ void scale (float *, int, int, float *);
__global__ void r1_update(float *, float *, int, int, int, int);
__host__ void gpuAssert(cudaError_t, char *, int);
__host__ void print_matrix(float *, int, int, int);
__host__ void print_help (int);
__host__ int option_parser(int, char **, int *, int *);

// Host code
int main (int argc, char **argv) {

    int m, n, lda, ldr, i;
    float *A_d, *R_d;
    //cudaEvent_t t_start, t_stop;
    // Get "m" and "n" from command line
    if (0 != option_parser(argc, argv, &m, &n)) {
        fprintf(stderr, "Can\'t continue, exiting now!\n"); 
        exit(EXIT_FAILURE);
    }
    size_t A_size = m * n * sizeof(float);
    size_t R_size = n * n * sizeof(float);

    lda = n; ldr = n;

    // Allocate input matrices A_h and R_h in host memory
    float *A_h = (float *) malloc(A_size);
    float *R_h = (float *) malloc(R_size);
    memset(R_h, 0, R_size);

    // Initialize input matrix
    for (i = 0; i < n; i++)
        A_h[i*lda + i] = i + 1;

    // Allocate matrices in device memory
    gpuErrchk (cudaMalloc(&A_d, A_size));
    gpuErrchk (cudaMalloc(&R_d, R_size));
    // Copy the A matrix from host memory to device memory
    gpuErrchk (cudaMemcpy(A_d, A_h, A_size, cudaMemcpyHostToDevice));
    // Set R matrix to 0 
    gpuErrchk (cudaMemset(R_d, 0, R_size));

    /**** Invoke kernel ****/
    dim3 dimBlock (THREAD_P_BLOCK, 1, 1);
    // dimGrid 'monodimensional' (just x value)
    dim3 dimGrid_M ((m + THREAD_P_BLOCK - 1) / THREAD_P_BLOCK, 1, 1);
    // dimGrid 'bidimensional' (x and y values)
    dim3 dimGrid_B (BLOCKS_P_GRID, (m + THREAD_P_BLOCK - 1) / THREAD_P_BLOCK,1);
    // Gram-Schmidt algorithm step by step
    for (i = 0; i < n; i++) {
        // Step #1 --> R(i,i:n-1) = A'(:,i) * A(:,i:n-1)
        xTA <<< dimBlock, dimGrid_B >>> (R_d, A_d, m, n, lda, ldr, i);
        // Step #3 (Is the scale of a column vector)
        scale <<< dimBlock, dimGrid_M >>> (A_d + i, m, lda, R_d + i*ldr + i);
        // Step #4 (Is the scale of a row)
        scale <<< dimBlock, dimGrid_M >>> (R_d + ldr*i, m, 1, R_d + i*ldr + i);
        // Step #5 --> A(:,i+1:n−1) = A(:,i+1:n−1) − A(:,i) ∗ R(i,i+1:n−1)
        r1_update <<< dimBlock, dimGrid_B >>> (A_d, R_d + i*lda, m, n, lda, i);
    }
    // Copy the results from device memory to host memory
    gpuErrchk (cudaMemcpy(A_h, A_d, A_size, cudaMemcpyDeviceToHost));
    // Free device memory
    cudaFree(A_d); cudaFree(R_d);
    // Free host memory
    free(A_h); free(R_h);

    return 0;
}

/**
 * ## Kernel 1
 * 
 * Rank 1 update of columns of A
 */
__global__ void r1_update (float *A, float *R, int m, int n, int lda, int k) {
    // get x,y cordinates
    unsigned int x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned int y = blockIdx.y * blockDim.y + threadIdx.x;

    if (x < m && y < n-k-1)
        A[x*lda + y + k + 1] -= A[x*lda + k] * R[y + k + 1];
}

/**
 * ## Kernel 2
 * 
 * matrix vector product
 * Performs R[i] =  x'A where x' is a row of A
 * 
 * How leading dimension is used for matrices: http://ibm.co/19PLtIX
 */
__global__ void xTA (float *R, float *A, int m, int n, int lda, int ldr, int k){
    // block column * block dim + column (computed by each thread)
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j;
    // upper triangular matrix 
    if (i < n - k) {
        for (j = 0; j < m; j++)
            R[k*ldr + k + i] += A[k*lda + j] * A[j*lda + k + i];
    }
}

/**
 * ## Kernel 3
 * 
 * mult. for constant s
 * d vector
 * ld leading dimension (distance from elements)
 */
__global__ void scale (float *d, int m, int ld, float *R_x) {
    // block colum * block dim + column (computed by each thread)
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    // s = sqrt(R(i,i))
    // Static initialization of shared variables is illegal in CUDA. 
    // The problem is that the semantics of how every thread should treat 
    // static initialization of shared memory is undefined in the 
    // programming model. Which thread should do the write? What happens if 
    // the value is not uniform between threads? How should the compiler 
    // emit code for such a case and how should the hardware run it?
    __shared__ float s; s = sqrt(*R_x);
    // and scale
    if (i < m) d[i*ld] /= s;
}

/*
 * GPU Error Handler (CUDA_SAFE_CALL deprecated from CUDA 5.0)
 */
__host__ void gpuAssert(cudaError_t code, char *file, int line) {

    if (code != cudaSuccess) {
        fprintf(stderr,"GPUassert: %s %s %d\n", 
            cudaGetErrorString(code), file, line);
        exit(code);
    }
}

/* 
 * Print matrix
 * 
 * Print a matrix passed as argument
 */
__host__ void print_matrix (float * matrix, int m, int n, int ld) {

    int i, j;

    for (i = 0; i < m; i++) {
        for (j = 0; j < n; j++)
            printf("%0.5f ", matrix[i*ld + j]);
        printf("\n");
    }
}

/*
 * The option parser
 *
 * The function parses the parameters passed from the command line and run 
 * their own procedures.
 *
 * Return value:
 *   0 on success
 *  -1 on failure
 *
 * Please, see http://www.gnu.org/software/libc/manual/html_node/Getopt.html
 * for further informations. (thanks to Frodo Looijaard)
 */
__host__ int option_parser (int argc, char **argv, int * m, int * n) {

    int opt;

    if (argc < 2) {
        fprintf(stderr, "The program needs arguments...\n\n");
        print_help(1);
    }

    opterr = 0;

    while ( -1 != (opt = getopt (argc, argv, "hr:c:"))) {
        switch (opt) {
            case 'h': 
                print_help(0);
            case 'r': 
                printf("optarg: %s\n", optarg); 
                if ((*m = atoi(optarg)) < 2) return -1; 
                break;
            case 'c': 
                printf("optarg: %s\n", optarg);
                if ((*n = atoi(optarg)) < 2 || *n > *m) return -1; 
                break;
            case '?':
                if (optopt == 'r' || optopt == 'c')
                    fprintf(stderr,"Option -%c requires an argument.\n",optopt);
                else if (isprint (optopt))
                    fprintf(stderr,"Unknown option `-%c'.\n", optopt);
                else
                    fprintf(stderr,"Unknown option chr `\\x%x'.\n", optopt);
                return -1;
            default:
                fprintf(stderr, "default switch-case statement reached\n");
                return -1;
        }
        //for (ii = optind; ii < argc; ii++)
        //  printf ("non-option argument %s\n", argv[ii]);
    }
    return 0;
}

/*
 * The helper
 *
 * Show the info to run the program in the correct way
 */
__host__ void print_help (int exit_code) {

    printf("\nPKG : %s\nPROGRAM : %s\nVERSION : %s\n\n",PACKAGE,PROG,VERSION);

    printf("%s [-h] [-r num of rows] [-c num of columns]\n\n", PROG);
    printf("  -h    print this help and exit\n");
    printf("  -r    provide the number of rows\n");
    printf("  -c    provide the number of colums\n\n");

    printf("  Example: ./qr_gpu -r 800 -c 600\n\n");

    exit_code == -1 ? exit(EXIT_FAILURE) : exit(EXIT_SUCCESS);
}

当我用 cuda-memcheck 运行程序时，我得到了这个结果:

[mcrociara@tesla project_CUDA]$ cuda-memcheck ./qr_gpu -r 4 -c 4 ========= CUDA-MEMCHECK optarg: 4 optarg: 4 GPUassert: unspecified launch failure src/qr_gpu.cu 99 ========= Invalid global read of size 4 ========= at 0x000000c8 in xTA ========= by thread (0,0,0) in block (0,0,0)

========= Address 0x3b5273104 is out of bounds

========= Invalid global read of size 4 ========= at 0x000000c8 in xTA ========= by thread (1,0,0) in block (0,0,0)

========= Address 0x3b5273108 is out of bounds

========= Invalid global read of size 4 ========= at 0x000000c8 in xTA ========= by thread (2,0,0) in block (0,0,0)

========= Address 0x3b527310c is out of bounds

========= ERROR SUMMARY: 3 errors

有人可能会帮助理解为什么？我在这个算法上实现了似乎可以正常工作的串行版本:

/*
 * QR decomposition via modified Gram-Schmidt algorithm
 */

// Libraries
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <ctype.h>
#include <time.h>
#include <math.h>

// Constant
#define PROG "QR_cpu"
#define VERSION "28.0913"
#define PACKAGE "QR-Decomposition"

// Prototypes
void r1_update(double *, double *, int);
void xTA (double *, double *, int);
void scale (double *, int, double);
void gram (double *, double *);
void print_matrix(double *);
int option_parser(int, char **);
void print_help (int);

// Number of rows
int M = -1; 
// Number of columns
int N = -1;
// Leading Dimension of A and R
int LDA = -1;
int LDR = -1;
// The clock
clock_t start, stop;

/*
 *
 */
int main (int argc, char **argv) {

    int i;

    // Get M, N from command line
    if (0 != option_parser(argc, argv)) {
        fprintf(stderr, "Bad option man!!!\n");
        exit(EXIT_FAILURE);
    }
    // Check the size of M and N 
    if (M > 5000 || N > 5000) {
        fprintf(stderr, "Too big man!!!\n");
        exit(EXIT_FAILURE);
    }

    // Set the leading dimension of A and R
    LDA = N; LDR = N;

    // Allocate memory for A and R 
    double *A = calloc(M * N, sizeof(*A));
    double *R = calloc(N * N, sizeof(*R));

    // Set the diagonal of A as A(i,i) = i + 1 with i from 0 to N−1
    for (i = 0; i < N; i++) 
        A[i*LDA + i] = i + 1;

    start = clock();
    gram(A, R);
    stop = clock();

    printf("\nTime: %0.4f\n\n",(stop - start) / (double)(CLOCKS_PER_SEC));
    // print_matrix(A);
    // print_matrix(R);

    free(A); free(R);

    return 0;
}

/**
 * Rank 1 update of columns of A
 */
void r1_update (double *A, double *R, int k) {

    int i, j;

    for(i = 0; i < M; i++)
        for(j = k + 1; j < N; j++)
            A[i*LDA + j] -= A[i*LDA + k] * R[j];
}

/**
 * Matrix vector product
 * Performs R[i] =  x'A where x' is a row of A
 * A : m x k, leading dimebsion, lda
 * 
 * How leading dimension is used for matrices: http://ibm.co/19PLtIX
 */
void xTA (double *R, double *A, int k) {

    int i, j;
    // upper triangular matrix 
    for (i = 0; i < N-k; i++)
        for (j = 0; j < M; j++)
            R[k*LDR + k + i] += A[k*LDA + j] * A[j*LDA + k + i];
}

/**
 * Mult. for constant s
 * d vector
 * ld leading dimension (distance from elements)
 */
void scale (double *d, int ld, double s) {

    int i;

    for (i = 0; i < M; i++) d[i*ld] *= s;
}

/**
 * Performs Modified Gram Schmidt
 * ortogonalization of columns of A
 * A m x n
 * R n x n
 */ 
void gram (double *A, double *R) {

    int i;
    double s;
    // Modified Gram Schmidt algorithm step by step
    for (i = 0; i < N; i++) {
        // Step #1 --> R(i,i:n-1) = A'(:,i) * A(:,i:n-1)
        xTA(R, A, i);
        // Step #2 (Normalizing) --> s = sqrt(R(i,i))
        s = 1 / sqrt(R[i*LDR + i]);
        // Step #3 (Is the scale of a column vector)
        scale(A + i, LDA, s);
        // Step #4 (Is the scale of a row)
        scale(R + LDR*i, 1, s);
        // Step #5 --> A(:,i+1:n−1) = A(:,i+1:n−1) − A(:,i) ∗ R(i,i+1:n−1)
        r1_update(A, R + i*LDA, i);
    }
}

/* 
 * Print Matrix
 * 
 * Print a matrix passed as argument
 */
void print_matrix (double * matrix) {

    int i, j;

    for (i = 0; i < M; i++) {
        for (j = 0; j < N; j++)
            printf("%0.2f ", matrix[i*LDA + j]);
        printf("\n");
    }
}

/*
 * The option parser
 *
 * The function parses the parameters passed from the command line and run 
 * their own procedures.
 *
 * Return value:
 *   0 on success
 *  -1 on failure
 *
 * Please, see http://www.gnu.org/software/libc/manual/html_node/Getopt.html
 * for further informations. (thanks to Frodo Looijaard)
 */
int option_parser (int argc, char **argv) {

    int opt;

    if (argc < 2) {
        fprintf(stderr, "This program needs arguments...\n\n");
        print_help(1);
    }

    opterr = 0;

    while ( -1 != (opt = getopt (argc, argv, "hr:c:"))) {
        switch (opt) {
            case 'h': 
                print_help(0);
            case 'r': 
                printf("optarg: %s\n", optarg); 
                if ((M = atoi(optarg)) < 2) return -1; 
                break;
            case 'c': 
                printf("optarg: %s\n", optarg);
                if ((N = atoi(optarg)) < 2 || N > M) return -1; 
                break;
            case '?':
                if (optopt == 'r' || optopt == 'c')
                    fprintf(stderr,"Option -%c requires an argument.\n",optopt);
                else if (isprint (optopt))
                    fprintf(stderr,"Unknown option `-%c'.\n", optopt);
                else
                    fprintf(stderr,"Unknown option chr `\\x%x'.\n", optopt);
                return -1;
            default:
                fprintf(stderr, "default switch-case statement reached\n");
                return -1;
        }
        //for (ii = optind; ii < argc; ii++)
        //  printf ("Non-option argument %s\n", argv[ii]);
    }
    return 0;
}

/*
 * The helper
 *
 * Shows the info to run the program in the correct way
 */
void print_help (int exit_val) {

    printf("\nPKG : %s\nPROGRAM : %s\nVERSION : %s\n\n",PACKAGE,PROG,VERSION);

    printf("%s [-h] [-r num of rows] [-c num of columns]\n\n", PROG);
    printf("  -h    print this help and exit\n");
    printf("  -r    provide the number of rows\n");
    printf("  -c    provide the number of colums\n\n");

    printf("  Example: ./qr_cpu -r 800 -c 600\n\n");

    exit_val == -1 ? exit(EXIT_FAILURE) : exit(EXIT_SUCCESS);
}

提前致谢!

文章推荐： c++ - 如何在 C++ 中逐字符输入？

文章推荐： node.js - Etherpad:连接身份验证的最佳方式

文章推荐： node.js - 使nodejs "require"路径独立？

文章推荐： c - for循环迭代问题

Vim - 如何使用 smartcase 进行/搜索，而使用 noic 进行 * 搜索？
我喜欢 smartcase，也喜欢 * 和 # 搜索命令。但我更希望 * 和 # 搜索命令区分大小写，而/和 ?搜索命令遵循 smartcase 启发式。是否有隐藏在某个地方我还没有找到的设置？我宁
通过 SSH 进行 SSH(或如何通过 SSH 进行 'proxify' SSH)
关闭。这个问题是off-topic .它目前不接受答案。想改进这个问题？ Update the question所以它是on-topic对于堆栈溢出。 10年前关闭。 Improve this qu
java - 使用一个 LDAP 进行 AD 身份验证失败，但使用另一 LDAP 进行 AD 身份验证通过
从以下网站，我找到了执行java AD身份验证的代码。 http://java2db.com/jndi-ldap-programming/solution-to-sslhandshakeexcepti
r - 在哪些情况下，人们更喜欢通过 reshape 进行 melt 而不是通过 plyr 进行 ddply？我正在努力学习它们，它们看起来很相似
似乎 melt 会使用 id 列和堆叠的测量变量 reshape 您的数据框，然后通过转换让您执行聚合。 ddply，从 plyr 包看起来非常相似..你给它一个数据框，几个用于分组的列变量和一个聚合
caching - 为什么 Facebook 在 memcached 中使用 TCP 进行 SET 和 UDP 进行 GET
我的问题是关于 memcached。 Facebook 使用 memcached 作为其结构化数据的缓存，以减少用户的延迟。他们在 Linux 上使用 UDP 优化了 memcached 的性能。 h
java - 在 Camel route ，使用 exec 组件使用 grep 进行 curl ，但使用 ${HOSTNAME} 进行 grep 无法正常工作
在 Camel route ，我正在使用 exec 组件通过 grep 进行 curl ，但使用 ${HOSTNAME} 的 grep 无法正常工作，下面是我的 Camel 路线。请在这方面寻求帮助。
django - 进行 "not in"查询
我正在尝试执行相当复杂的查询，在其中我可以排除与特定条件集匹配的项目。这是一个 super 简化的模型来解释我的困境: class Thing(models.Model) user = mod
django - 进行 "not in"查询
我正在尝试执行相当复杂的查询，我可以在其中排除符合特定条件集的项目。这里有一个 super 简化的模型来解释我的困境: class Thing(models.Model) user = mod
angular - 进行 Angular 内容投影的现代方法是什么？
我发现了很多嵌入/内容项目的旧方法，并且我遵循了在这里找到的最新方法(我假设):https://blog.angular-university.io/angular-ng-content/ 我正在尝试
使用 NextJS 进行 Fastify
我正在寻找如何使用 fastify-nextjs 启动 fastify-cli 的建议我曾尝试将代码简单地添加到建议的位置，但它不起作用。 'use strict' const path = req
javascript - 进行 gatsby 构建时未定义窗口
我正在尝试将振幅 js 与 React 和 Gatsby 集成。做 gatsby developer 时一切看起来都不错，因为它发生在浏览器中，但是当我尝试 gatsby build 时，我收到以下错
java - 进行 Null 检查的频率和位置
我试图避免过度执行空值检查，但同时我想在需要使代码健壮的时候进行空值检查。但有时我觉得它开始变得如此防御，因为我没有实现 API。然后我避免了一些空检查，但是当我开始单元测试时，它开始总是等待运行时异
使用 NOT 进行 Kibana 搜索
尝试进行包含一些 NOT 的 Kibana 搜索，但获得包含 NOT 的结果，因此猜测我的语法不正确: "chocolate" AND "milk" AND NOT "cow" AND NOT "tr
iphone - 进行 Facebook 集成时出错
我正在使用开源代码共享包在 iOS 中进行 facebook 集成，但收到错误“FT_Load_Glyph failed: glyph 65535: error 6”。我在另一台 mac 机器上尝试了
r - 进行 Tobit 回归时的奇异性错误
我正在尝试估计一个标准的 tobit 模型，该模型被审查为零。变量是因变量 : 幸福自变量 : 城市(芝加哥，纽约)，性别(男，女)，就业(0=失业，1=就业)，工作类型(失业，蓝色，白色
从多个文件夹中对多个 jar 进行 gradle
我有一个像这样的项目布局样本/ 一种/ 源/ 主要的/ java / java 资源/ .jpg 乙/ 源/ 主要的/ java / B.java 资源/ B.jpg 构建.gradle 设置.gr
javascript - 进行 fetch 调用时出错
如何循环遍历数组中的多个属性以及如何使用map函数将数组中的多个属性显示到网页 import React, { Component } from 'react'; import './App.css'
javascript - 进行 AJAX 调用时加载程序不显示
我有一个 JavaScript 函数，它进行 AJAX 调用以返回一些数据，该调用是在选择列表更改事件上触发的。我尝试了多种方法来在等待时显示加载程序，因为它当前暂停了选择列表，从客户的 Angul
java - 进行 null 检查的更短方法
可能以前问过，但找不到。我正在用以下形式写很多语句: if (bar.getFoo() != null) { this.foo = bar.getFoo(); } 我想到了三元运算符，但我认
javascript - 进行 JavaScript 验证后短信消失
我有一个表单，在将其发送到 PHP 之前我正在执行一些验证 JavaScript，验证后的 JavaScript 函数会发布用户在中输入的文本。页面底部的标签；然而，此消息显示短暂，然后消失...

太空宇宙

个人简介

我是一名优秀的程序员,十分优秀！

作者热门文章

滴滴打车优惠券免费领取

全站热门文章

首页

博学

6Ren·AI

商城

c - 如何在 C 和 CUDA 中通过修改后的 Gram-Schmidt 方法进行 QR 分解

========= Address 0x3b5273104 is out of bounds

========= Address 0x3b5273108 is out of bounds

========= Address 0x3b527310c is out of bounds