Cuda: resolución de mínimos cuadrados, baja velocidad

Recientemente, uso Cuda para escribir un algoritmo llamado 'búsqueda de correspondencia ortogonal'. En mi código feo de Cuda, la iteración completa tarda 60 segundos, y Eigen lib tarda solo 3 segundos ...

En mi código, la matriz A es [640,1024] e y es [640,1], en cada paso selecciono algunos vectores de A para componer una nueva matriz llamada A_temp [640, itera], iter = 1: 500. Nueva una matriz MaxDex_Host [] en la CPU para decir qué columna seleccionar.

Quiero obtener x_temp [itera, 1] de A_temp * x_temp = y usando el mínimo cuadrado, uso una API de cula 'culaDeviceSgels' y una API de multiplicación de vector de matriz de cublas.

Entonces culaDeviceSgels llamaría 500 veces, y creo que esto sería más rápido que el QR.Sovler de Eigen lib.

Compruebo el análisis de rendimiento de Nisight, descubrí que la historia del sueño lleva mucho tiempo. Inicializo cublas antes de la iteración y lo desecho después de obtener el resultado. Entonces, quiero saber qué es el Custreamdestory, diferente con Cublasdestory.

El principal problema es memcpy y la función 'gemm_kernel1x1val'. Creo que esta función es de 'culaDeviceSgels'

while (itera <500): uso cublasSgemv y cublasIsamax para obtener MaxDex_Host [itera], luego

        MaxDex_Host[itera]=pos;
    itera++; 
    float* A_temp_cpu=new float[M*itera]; // matrix all in col-major
    for (int j=0;j<itera;j++) // to  get A_temp [M,itera] , the MaxDex_Host[] shows the positon of which column of A to chose , 
    {
        for (int i=0;i<M;i++) //M=640 , and A is 640*1024 ,itera is add 1 each step
        {
            A_temp_cpu[j*M+i]=A[MaxDex_Host[j]*M+i];
        }
    }
          // I must allocate one more array because culaDeviceSgels will decompose the one input Array ,  and I want to use A_temp after least-square solving.
    float* A_temp_gpu;
    float* A_temp2_gpu;  
    cudaMalloc((void**)&A_temp_gpu,Size_float*M*itera);
    cudaMalloc((void**)&A_temp2_gpu,Size_float*M*itera);
    cudaMemcpy(A_temp_gpu,A_temp_cpu,Size_float*M*itera,cudaMemcpyHostToDevice);
    cudaMemcpy(A_temp2_gpu,A_temp_gpu,Size_float*M*itera,cudaMemcpyDeviceToDevice);
    culaDeviceSgels('N',M,itera,1,A_temp_gpu,M,y_Gpu_temp,M);// the x_temp I want is in y_Gpu_temp's return value ,  stored in the y_Gpu_temp[0]——y_Gpu_temp[itera-1]
     float* x_temp;
    cudaMalloc((void**)&x_temp,Size_float*itera);
    cudaMemcpy(x_temp,y_Gpu_temp,Size_float*itera,cudaMemcpyDeviceToDevice);

El manejo de la memoria de Cuda parece demasiado complejo, ¿hay algún otro método conveniente para resolver el mínimo cuadrado?

Respuestas a la pregunta(1)

Su respuesta a la pregunta