warp 上的 Jacobi 旋转法计算 Hermitian 矩阵特征值和特征向量的 cuda 实现

于 2024-04-11 12:27:16 发布
阅读量224 收藏
点赞数 1
文章标签： 矩阵 c++ 算法 cuda
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/eloudy/article/details/137634253
版权
本文详细介绍了在CUDA平台上使用cuComplex数据类型实现Jacobi旋转算法，包括Schur变换、右/左旋转操作以及Hermitian矩阵的处理，同时展示了如何调用cuSolver进行Hermitianeigendecomposition的实际应用。
摘要由CSDN通过智能技术生成

1，原理

Jacobi旋转的 cuComplex数据类型上的实现

2，实现代码


#include <cuda_runtime.h>
#include <cusolverDn.h>
#include <cuComplex.h>
#include <iostream>
#include <string>


#define ZERO_IJ 1
using namespace std;

__host__ __device__
float2 zero__ = {0.0f, 0.0f};

__host__ __device__
float2  zero_f2 = {0.0f, 0.0f};



#if PRINT_DEBUG

template <typename T1>
void print_vector(T1 *A, int N)
{
        for (int idx = 0; idx < N; idx++)
        {
                // printf("%f\n", A[idx]);
                cout << A[idx] << endl;
        }
        cout << endl;
}

void print_Int_vector(int *A, int N)
{
        for (int idx = 0; idx < N; idx++)
        {
                printf("%5d ", A[idx]);
        }
}

template <typename T1>
void print_matrix(T1 *A, int M, int N, int lda)
{
        for (int i = 0; i < M; i++)
        {
                for (int j = 0; j < N; j++)
                {
                        // printf("%7.4f ", A[i + j*lda]);
                        cout << A[i + j * lda] << " ";
                }
                // printf("\n");
                cout << endl;
        }
}
template <typename T2>
void print_matrix_complex(T2 *A, int M, int N, int lda)
{
        for (int i = 0; i < M; i++)
        {
                for (int j = 0; j < N; j++)
                { //*(sizeof(complex<float>)/sizeof(float))
                        T2 tmp;
                        tmp = A[i + j * lda];
                        cout << "(" << tmp.x << " + " << tmp.y << "*j) "
                        // printf("(%7.4f, %7.4f)", tmp.x, tmp.y);
                }
                cout << endl;
                // printf("\n");
        }
}

#endif

/

__device__ float signf_thread(float a)
{
        union f_i
        {
                float a;
                int x;
        };
        union f_i tmp;
        tmp.a = a;

        return (tmp.x >> 31) ? -1.0f : 1.0f;
}

__device__ void Schur2_Hermite_thread(float Aii, float Ajj, float2 Aij, float norm_Aij, float2 &cC, float2 &s_plus, float2 &s_neg, float &t)
{
        float Aii_Ajj = Aii - Ajj;
        t = 2.0f * norm_Aij * signf_thread(Aii_Ajj) / (fabs(Aii_Ajj) + sqrt(Aii_Ajj * Aii_Ajj + 4.0f * norm_Aij * norm_Aij));
        float c = 1.0f / sqrt(1 + t * t);
        float s = c * t;
        float sDNAij; // s/|Aij|

        sDNAij = s / norm_Aij;
        s_plus = make_cuFloatComplex(sDNAij * Aij.x, sDNAij * Aij.y);
        s_neg = cuConjf(s_plus);
        cC = make_cuFloatComplex(c, 0.0f);
}

__device__ void right_rotate_H_thread(float2 *A, int n, int lda, float2 cC, float2 s_neg, float2 s_plus, int i, int j)
{
        // for (int r = 0; r < n; r++) {// right_rotate  i-th column  j-th column  changed
        unsigned r = threadIdx.x;
        float2 Ari, Arj;

        Ari = A[r + i * lda];
        Arj = A[r + j * lda];
        // if(r!=i)
        A[r + i * lda] = cuCaddf(cuCmulf(Ari, cC), cuCmulf(Arj, s_neg)); // Ari * c    + Arj * S-
        // if(r!=j)
        A[r + j * lda] = cuCsubf(cuCmulf(Arj, cC), cuCmulf(Ari, s_plus)); // Ari *(-S+) + Arj * c
        //}
}

__device__ void left_rotate_H_thread(float2 *A, int n, int lda, float2 cC, float2 s_neg, float2 s_plus, int i, int j)
{
        // for (int r = 0; r < n; r++) {// left_rotate i-th row   j-th row changed
        unsigned r = threadIdx.x;
        float2 Air, Ajr;

        Air = A[i + r * lda];
        Ajr = A[j + r * lda];

        A[i + r * lda] = cuCaddf(cuCmulf(cC, Air), cuCmulf(s_plus, Ajr)); //  c * Air   + S+ * Ajr
        A[j + r * lda] = cuCsubf(cuCmulf(cC, Ajr), cuCmulf(s_neg, Air));  // -S-* Air   +  c * Ajr
        //}
}

__device__ void set_V_2_E_thread(unsigned n, float2 *V_d, int ldv)
{
        float2 one_ = make_cuFloatComplex(1.0f, 0.0f);
        float2 zero_ = make_cuFloatComplex(0.0f, 0.0f);

        unsigned tid = threadIdx.x;

        for (unsigned j = 0; j < n; j++)
        {
                V_d[tid + j * ldv] = zero_;
        }

        V_d[tid + tid * ldv] = one_;
}

__device__ void sort_thread(float *sig, int N)
{

        int x = threadIdx.x;
        extern __shared__ float s_a[];
        if (x < N)
                s_a[x] = sig[x];
        __syncthreads();

        for (int i = 0; i < N; i++)
        {

                int j = i % 2;
                int idx = 2 * x + j;
                if (idx + 1 < N && s_a[idx] > s_a[idx + 1])
                {
                        float tmp = s_a[idx];
                        s_a[idx] = s_a[idx + 1];
                        s_a[idx + 1] = tmp;
                }
                __syncthreads();
        }
        if (x < N)
                sig[x] = s_a[x];
        __syncthreads();
}

__device__ void sort_value_id_thread(float *sig, int N, float *s_a, unsigned *ids)
{

        int x = threadIdx.x;
        // extern __shared__ float s_a[];
        if (x < N)
        {
                s_a[x] = sig[x];
                ids[x] = x;
        }
        __syncthreads();

        for (int i = 0; i < N; i++)
        {

                int j = i % 2;
                int idx = 2 * x + j;
                if (idx + 1 < N && s_a[idx] > s_a[idx + 1])
                {
                        float tmp = s_a[idx];
                        s_a[idx] = s_a[idx + 1];
                        s_a[idx + 1] = tmp;
                        unsigned tmpid = ids[idx];
                        ids[idx] = ids[idx + 1];
                        ids[idx + 1] = tmpid;
                }
                __syncthreads();
        }
        if (x < N)
                sig[x] = s_a[x];
        __syncthreads();
}

__global__ void hermite_jacobi_kernel(float2 *A, int n, int lda, bool cal_Vectors, float2 *V, int ldv, float *Sig, float eps, unsigned jt)
{
        unsigned tid = threadIdx.x;
        if (tid < n)
                set_V_2_E_thread(n, V, ldv);

        unsigned count = 0;

        while (count < jt)
        {
                count++;
                for (int i = 0; i < n - 1; i++)
                {
                        for (int j = i + 1; j < n; j++)
                        {

                                float s, c, t;
                                float2 Aij;
                                float norm_Aij;
                                float Aii, Ajj;
                                float2 s_plus, s_neg;
                                float2 cC;

                                Aii = A[i + i * lda].x;
                                Ajj = A[j + j * lda].x;
                                Aij = A[i + j * lda];

                                norm_Aij = cuCabsf(Aij);
                                if (norm_Aij < eps)
                                        continue;

                                Schur2_Hermite_thread(Aii, Ajj, Aij, norm_Aij, cC, s_plus, s_neg, t);

                                if (tid < n)
                                {
                                        right_rotate_H_thread(A, n, lda, cC, s_neg, s_plus, i, j);
                                        left_rotate_H_thread(A, n, lda, cC, s_neg, s_plus, i, j);
                                }
#if ZERO_IJ
                                if (tid == 63)
                                { // third warp
                                        A[i + j * lda] = zero__;
                                        A[j + i * lda] = zero__;
                                        // float A_i_i = c * c * Aii + 2.0f * c * s * norm_Aij + s * s * Ajj;
                                        //  printf("A[i,i] %d= %f\n", count, A[i+i*lda].x);
                                        //  printf("A[i,i] %d= %f\n", count, A_i_i);
                                        float a_ii_look = Aii + t * norm_Aij;
                                        // printf("a_ii_ %d = %f\n", count, a_ii_look);

                                        // float A_j_j = s * s * Aii - 2.0f * c * s * norm_Aij + c * c * Ajj;
                                        //  printf("A[j,j] %d= %f\n", count, A[j + j * lda].x);
                                        //  printf("A[j,j] %d= %f\n", count, A_j_j);
                                        float a_jj_look = Ajj - t * norm_Aij;
                                        // printf("a_jj_ %d = %f\n", count, a_jj_look);

                                        // A[i + i * lda].x = A_i_i;        A[i + i * lda].y = 0.0f;
                                        // A[j + j * lda].x = A_j_j;        A[j + j * lda].y = 0.0f;
                                        A[i + i * lda].x = a_ii_look;
                                        A[i + i * lda].y = 0.0f;
                                        A[j + j * lda].x = a_jj_look;
                                        A[j + j * lda].y = 0.0f;
                                }
                //__syncthreads();
//                if(threadIdx.x == 0) printf("count = %d\n", count);
#endif
                                if (cal_Vectors)
                                { // another warp do this
                                        if (tid < n)
                                                right_rotate_H_thread(V, n, ldv, cC, s_neg, s_plus, i, j);
                                }
                //__syncthreads();
                        }

                }
        }

        // for (int dia = 0; dia < n; dia++) {
        if (tid < n)
        {
                Sig[tid] = A[tid + tid * lda].x;
        }

#if 1
        // sort Sigma
        int N = n;
        unsigned x = tid; // threadIdx.x;
        __shared__ unsigned ids[64];
        __shared__ float s_a[64];

        sort_value_id_thread(Sig, N, s_a, ids);

        for (int col = 0; col < n; col++)
        {
                if (tid < n)
                {
                        A[tid + col * lda] = V[tid + ids[col] * ldv];
                }

        }
#endif
}

void gpu_hermite_jacobi(float2 *A_d, int n, int lda, bool cal_Vectors, float2 *V_d, int ldv, float *Sig, float eps, unsigned jt)
{
        dim3 grid_;
        dim3 block_;

        block_.x = 32;
        //__global__ void hermite_jacobi_kernel(float2* A, int n, int lda, bool cal_Vectors, float2* V, int ldv, float* Sig, float eps, unsigned jt)
#if MEASURE_TIME
        cudaEvent_t start, end;
        cudaEventCreate(&start);
        cudaEventCreate(&end);
        //
        cudaEventRecord(start);
#endif
        hermite_jacobi_kernel<<<grid_, block_, 0, NULL>>>(A_d, n, lda, cal_Vectors, V_d, ldv, Sig, eps, jt);

#if MASURE_TIME
        cudaEventRecord(end);
        cudaEventSynchronize(end);
        //
        float time_ms = 0.0f;
        cudaEventElapsedTime(&time_ms, start, end);
        std::cout << "CUDA Kernel time: " << time_ms << " ms" << std::endl;

        cudaEventDestroy(start);
        cudaEventDestroy(end);
        std::cout << "end gpu 1" << std::endl;
#endif
}

cusolverStatus_t cusolverDnCheevj_(
        cusolverEigMode_t jobz, // CUSOLVER_EIG_MODE_NOVECTOR     CUSOLVER_EIG_MODE_VECTOR
        cublasFillMode_t uplo,  // CUBLAS_FILL_MODE_LOWER    CUBLAS_FILL_MODE_UPPER     CUBLAS_FILL_MODE_FULL
        int N,                                  // n
        cuComplex *A,                   // A
        int lda,                                // lda
        float *W,                               // float *W= nullptr;                   cudaMalloc((void**)&W, N*sizeof(float));
        cuComplex *work,                // cuComplex *work=nullptr;             cudaMalloc((void**)&work, lwork);
        int lwork,                              // 3*N*(N+5)
        int *info)                              // 0: fine;   -i: the i-th parameter failed;   i:  i-th off-diagonal element cannot converge to zero;
{
        // return cusolverDnXXevd<float, cuComplex>(handle, jobz, uplo, N, A, lda, W, work, lwork, info);
        cudaMemset(info, 0x00, sizeof(int));
        bool cal_Vectors = (jobz == CUSOLVER_EIG_MODE_VECTOR);
        int ldv = N;
    //void gpu_hermite_jacobi(float2 *A_d, int n, int lda, bool cal_Vectors, float2 *V_d, int ldv, float *Sig, float eps, unsigned jt)
        gpu_hermite_jacobi(A, N, lda, cal_Vectors, work, ldv, W, 0.00000001f, 7);
        //new gpu_Cheevd_Jacobi(N, A, lda, cal_Vectors, work, ldv, W, 0.00001f,70);

        return CUSOLVER_STATUS_SUCCESS;
}









// CUDA API error checking
#define CUDA_CHECK(err)                                                                            \
    do {                                                                                           \
        cudaError_t err_ = (err);                                                                  \
        if (err_ != cudaSuccess) {                                                                 \
            printf("CUDA error %d at %s:%d\n", err_, __FILE__, __LINE__);                          \
            throw std::runtime_error("CUDA error");                                                \
        }                                                                                          \
    } while (0)

// cusolver API error checking
#define CUSOLVER_CHECK(err)                                                                        \
    do {                                                                                           \
        cusolverStatus_t err_ = (err);                                                             \
        if (err_ != CUSOLVER_STATUS_SUCCESS) {                                                     \
            printf("cusolver error %d at %s:%d\n", err_, __FILE__, __LINE__);                      \
            throw std::runtime_error("cusolver error");                                            \
        }                                                                                          \
    } while (0)



void print_complex_matrix(cuComplex* A, int m, int n, int lda)
{
    for(int i=0; i<(m<12? m:12); i++){
        for(int j=0; j<(n<12? n:12); j++){
            printf("(%6.3f+(", A[i + j*lda].x);
            printf("%6.3f*j)), ", A[i + j*lda].y);
        }
        printf("; ...\n");
    }
}

void init_Hermitian_matrix(cuComplex *A, int m, int n ,int lda, int seed)
{
    srand(seed);

    for(int i=0; i<lda; i++){
        for(int j=0; j<m; j++){
            if(i<=j){
                A[i + j*lda].x = (rand()%2000)/1000.0f;
                A[i + j*lda].y = (rand()%2000)/1000.0f;
                if(i==j)
                    A[i + j*lda].y = 0.0f;
            }else{
                A[i + j*lda].x = A[j + i*lda].x;
                A[i + j*lda].y = -A[j + i*lda].y;
            }
        }
    }
}

void complex_gemm_NT(cuComplex *A, int lda, cuComplex *B, int ldb, cuComplex *C, int ldc, int M, int N, int K)
{
    cuComplex zero_c;
    zero_c.x = 0.0f;
    zero_c.y = 0.0f;

    for(int i=0; i<M; i++){
        for(int j=0; j<N; j++){
            cuComplex sigma = zero_c;
            for(int k=0; k<K; k++){
                sigma = cuCaddf(sigma, cuCmulf(A[i + k*lda], cuConjf(B[j+ k*ldb])));
            }
            C[i + j*ldc] = sigma;
        }
    }
}










int main(int argc, char *argv[]) {
    cusolverDnHandle_t cusolverH = NULL;
    cudaStream_t stream = NULL;
    syevjInfo_t syevj_params = NULL;

    const int m = 7;
    const int lda = m;

    cuComplex *A = nullptr;
    A = (cuComplex*)malloc(lda*m*sizeof(cuComplex));

    cuComplex *V = nullptr;
    V = (cuComplex*)malloc(lda*m*sizeof(cuComplex));
    float *W = nullptr;
    W = (float*)malloc(m*sizeof(float));

    init_Hermitian_matrix(A, m, m, lda, 2024);

    cuComplex *d_A = nullptr;
    float *d_W = nullptr;
    int *devInfo = nullptr;
    cuComplex *d_work = nullptr;
    int lwork = 4096;
    int info_gpu = 0;

    /* configuration of syevj  */
    const double tol = 1.e-7;
    const int max_sweeps = 15;
    const cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvectors.
    //const cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
    const cublasFillMode_t uplo = CUBLAS_FILL_MODE_UPPER;

    /* numerical results of syevj  */
    double residual = 0;
    int executed_sweeps = 0;

    printf("tol = %E, default value is machine zero \n", tol);
    printf("max. sweeps = %d, default value is 100\n", max_sweeps);

    printf("A = (matlab base-1)\n");
    //print_matrix(m, m, A, lda);
    print_complex_matrix(A, m, m, lda);
    printf("=====\n");

    /* step 3: copy A to device */
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_A), sizeof(cuComplex) * lda * m));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_W), sizeof(float) * m));
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&devInfo), sizeof(int)));

    CUDA_CHECK(cudaMemcpy(d_A, A, sizeof(cuComplex) * lda * m, cudaMemcpyHostToDevice));

    printf("LL:: lwork = %d\n", lwork);
    CUDA_CHECK(cudaMalloc(reinterpret_cast<void **>(&d_work), sizeof(cuComplex) * lwork));

    /* step 5: compute eigen-pair   */
    CUSOLVER_CHECK(cusolverDnCheevj_( jobz, uplo, m, d_A, lda, d_W, d_work, lwork, devInfo));
    cudaDeviceSynchronize();
    CUDA_CHECK(cudaMemcpy(V, d_A, sizeof(cuComplex) * lda * m, cudaMemcpyDeviceToHost));
    CUDA_CHECK(cudaMemcpy(W, d_W, sizeof(float) * m, cudaMemcpyDeviceToHost));
    CUDA_CHECK(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

    //CUDA_CHECK(cudaStreamSynchronize(stream));

    if (0 == info_gpu) {
        printf("syevj converges \n");
    } else if (0 > info_gpu) {
        printf("%d-th parameter is wrong \n", -info_gpu);
        exit(1);
    } else {
        printf("WARNING: info = %d : syevj does not converge \n", info_gpu);
    }

    printf("Eigenvalue = (matlab base-1), ascending order\n");
    for (int i = 0; i < m; i++) {
        printf("W[%d] = %E\n", i + 1, W[i]);
    }
#if 1
    printf("V = (matlab base-1)\n");
    print_complex_matrix(V, m, m, lda);
    printf("=====\n");
#endif

    cuComplex *E = nullptr;
    E = (cuComplex*)malloc(m*m*sizeof(cuComplex));

    //void complex_gemm_NT(cuComplex *A, int lda, cuComplex *B, int ldb, cuComplex *C, int ldc, int M, int N, int K)
    complex_gemm_NT(V, lda, V, lda, E, m, m, m, m);
    printf("E =\n");
    print_complex_matrix(E, m, m, m);

    printf("residual |A - V*W*V**H|_F = %E \n", residual);
    printf("number of executed sweeps = %d \n", executed_sweeps);

    /* free resources */
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_W));
    CUDA_CHECK(cudaFree(devInfo));
    CUDA_CHECK(cudaFree(d_work));
    CUDA_CHECK(cudaDeviceReset());

    return EXIT_SUCCESS;
}

3， 编译

nvcc hello_jacobi_hermitian.cu

nvcc hello_jacobi_hermitian.cu -L /usr/local/cuda/lib64 -lcudart

4， 运行

./a.out

Eloudy
关注
  • 1
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
CUDA：线性系统a*x=b上的Jacobi方法实例
希望我的博客，能帮上你解决学习中工作中所遇到的问题
12-21 151
CUDA：线性系统a*x=b上的Jacobi方法实例
Jacobi旋转法
04-26
是真的jacobi旋转法 c语言的 。不骗人。
雅克比旋转
开心一刻的专栏
11-13 2405
在数值线性代数中，雅可比旋转是 n 维内积空间的二维线性子空间的旋转 Qkℓ，在用做相似变换的时候，被选择来置零 n×n 实数对称矩阵 A 的非对角元素的对称对: 它是雅可比特征值算法的核心运算，它是数值上稳定的并适合用并行计算实现。 注意到只有 A 的行 k 和 ℓ 与列 k 和 ℓ 受到影响，并且 A′ 将保持对称。还有给 Qkℓ 的明显的矩阵很少被计算，转而计算辅助
雅可比旋转（Jacobi法）求对称矩阵的特征值和特征向量
妖米的博客
08-24 7245
该方法是求解对称矩阵全部特征值和特征向量的一种方法，它基于以下结论： ①任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q,使得 QTAQ=diag(λ1,λ2,…,λn)Q^TAQ=diag(λ1,λ2,…,λn)QTAQ=diag(λ1,λ2,…,λn)②在正交相似变换下,矩阵元素的平方和不变。 即设A=(aij)n×nA=(a_{ij})_{n×n}A=(aij​)n×n​,Q为正交矩阵,记B=QTAQ=(bij)n×nB=Q^TAQ=(b_{ij})_{n×n}B=QTAQ=(bij​)
【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（三）：Jacobi 旋转法【理论到程序】
天地玄黄 魑魅魍魉 风花雪月 商角徵羽 暂时停更 十月重见
11-30 1845
本文介绍了矩阵特征值与特征向量的计算 之Jacobi 旋转法的基本原理和步骤，通过具体的矩阵示例演示其应用过程，并给出其Python实现。
【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（二）：Jacobi 过关法（Jacobi 旋转法的改进）【理论到程序】
天地玄黄 魑魅魍魉 风花雪月 商角徵羽 暂时停更 十月重见
12-01 1662
本文介绍了矩阵特征值与特征向量的计算之Jacobi 过关法的基本原理和步骤，并给出其Python实现。
GPU上大规模稀疏矩阵特征值计算高效算法之一——GPU介绍
bob908的专栏
08-19 6766
1.GPU基本知识 2.稀疏矩阵存储方式 3.ju z
cuda矩阵相乘_CUDA计算矩阵相乘
weixin_39793319的博客
12-20 598
1.最简单的 kernel 函数__global__ void MatrixMulKernel( float* Md, float* Nd, float* Pd, int Width){int tx = threadIdx.x; // cloumnint ty = threadIdx.y; // rowfloat Pvalue = 0;for (int k = 0; k{float Mdele =...
基于CUDA的GPU计算PI值【100011443】
04-03
在项目"basedon_cuda_gpu_calculate_pi"中，我们可以期待看到如何实现CUDA程序来并行计算π值，包括初始化设备、分配内存、定义计算函数、同步线程、读写结果等步骤。此外，可能会有不同线程分块配置的实验，以及对...
warp-rnnt:CUDA-Warp RNN-换能器
05-06
CUDA-Warp RNN-换能器 一个GPU实现RNN换能器（格雷夫斯 ， ）。 该代码是从移植的（由Awni Hannun编写），并且充分利用了CUDA warp机制。 损失的主要瓶颈是基于动态规划算法的前进/后退。 特别是，存在一个嵌套循环...
jacobi旋转法的VB实现
Candy__1的博客
01-22 973
jacobi旋转法的VB实现一、jacobi方法概要1.jacobi方法的基本思想:jacobi旋转法是通过一组平面旋转变换（构建一个平面旋转矩阵）对**实对称方阵A**进行相似对角化，化其为*对角阵*，进而求出特征值与特征向量的方法。2.jacobi方法的基本步骤<1>.选取非对角线元素的主元素(这里只画出了上三角阵）ai(p,q)a_i(p,q)ai​(p,q)<2> 构造如下形式的旋转矩阵其中：(1) app=aqqa_{pp}=a_{qq}app​=aqq​时，tan⁡(2θ
配准中雅可比行列式（Jacobian determinant）的原理、意义、代码实现以及可视化（5.14更新）
Timmymm的博客
02-19 1万+
2019-5-14更新 最近两天我研究了一下雅可比行列式的问题，重新检查了别人的公开代码与可视化，发现了一些问题，现总结如下： 检查雅可比行列式的计算原理与其代码实现，我发现之前的代码存在问题，现在进行更新纠正； 关于其代码实现，我还有一个疑问没有解答，这可能与其数学原理有关，需要进一步的验证； 对雅可比行列式的计算结果进行可视化，本文提供了代码参考，并进行了解释，但还没解决。 1.正确的代...
【数值计算方法（黄明游）】矩阵特征值与特征向量的计算（五）：Householder方法【理论到程序】
最新发布
天地玄黄 魑魅魍魉 风花雪月 商角徵羽 暂时停更 十月重见
12-02 1706
本文介绍了矩阵特征值与特征向量的计算之Householder方法的基本原理和步骤，包括例题解析及Python实现。
光流法和直接法介绍
m0_37604894的博客
10-09 851
光流为追踪源图像某个点在其他图像中的运动， 一般分为稀疏光流和稠密光流，其本质是估计像素在不同时刻图像中的运动本文主要介绍LK光流。
Jacobi 方法
guanglihuan
11-05 2238
第三节 Jacobi 方法     Jacobi方法是求对称矩阵的全部特征值以及相应的特征向量的一种方法，它是基于以下两个结论     1) 任何实对称矩阵A可以通过正交相似变换成对角型，即存在正交矩阵Q,使得                 QT AQ = diag(λ1 ,λ2 ,…,λn ) (3.1) 其中λi(i=1,2,…,n)是A的特征值，Q中各列为相应的特征向量
基于单边jacobi的奇异值分解（SVD）
k531623594的专栏
02-03 1万+
基于单边jacobi的奇异值分解（SVD）对于奇异值分解（SVD），相信很多学过线性代数/高等代数的同学都不会很陌生，但是怎么实现呢？接下来就来详细说说。　　为了方便讨论，本文所有的讨论仅限定于实数空间。奇异值分解的含义就是将矩阵A分解成一个酉矩阵U，一个准对角矩阵S还有一个酉矩阵V。首先要说明一点，SVD分解是存在的但不唯一，这个有兴趣的读者可以思考一下。利用数学公式可以写成：  A=USV′
SLAM代码（光流法）
热门推荐
何必浓墨重彩
09-11 1万+
光流法的目标是完成图像点的跟踪, 因此这里假设存在一个输入图像I, 以及要跟踪的点x, 存在另外一个图像块T, 我们的目标是完成图像块T到输入图像I的匹配. 文章针对Lucas-Canade光流法做了一个总结, 文章对lucas-canade (Forward Additive, FA)算法做了简介, 引入了Compositional算法以及Inverse方法. 因此对应组合形成4种方法分别是
并行计算奇异值分解--Jacobi旋转(转)
qq_37065641的博客
10-24 1200
https://wenku.baidu.com/view/38bab6c90408763231126edb6f1aff00bed57045.html http://blog.sina.com.cn/s/blog_4ff49c7e0102vysp.html
并行计算奇异值分解--Jacobi旋转
weixin_34318956的博客
07-06 900
鉴于矩阵的奇异值分解SVD在工程领域的广泛应用（如数据压缩、噪声去除、数值分析等等，包括在NLP领域的潜在语义索引LSI核心操作也是SVD），今天就详细介绍一种SVD的实现方法--Jacobi旋转法。跟其它SVD算法相比，Jacobi法精度高，虽然速度慢，但容易并行实现。 一些链接  http://cdmd.cnki.com.cn/Article/CDMD-10285-1012286387.h...
GPU-Warp驱动的彩色法有限元矩阵生成与高效装配
本文主要探讨了"基于GPU-Warp的着色法有限元矩阵生成与装配"这一主题，它发表在《计算设计与工程》杂志的2019年第六期，卷号705。作者Utpal Kiran、Deepak Sharma Chang和Sachin Singh Gautam来自印度理工学院机械...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值