/* Copyright (c) 2022, NVIDIA CORPORATION. All rights reserved. * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * * Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in the * documentation and/or other materials provided with the distribution. * * Neither the name of NVIDIA CORPORATION nor the names of its * contributors may be used to endorse or promote products derived * from this software without specific prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY * OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ /* * This sample implements a preconditioned conjugate gradient solver on * the GPU using CUBLAS and CUSPARSE. Relative to the conjugateGradient * SDK example, this demonstrates the use of cusparseScsrilu02() for * computing the incompute-LU preconditioner and cusparseScsrsv2_solve() * for solving triangular systems. Specifically, the preconditioned * conjugate gradient method with an incomplete LU preconditioner is * used to solve the Laplacian operator in 2D on a uniform mesh. * * Note that the code in this example and the specific matrices used here * were chosen to demonstrate the use of the CUSPARSE library as simply * and as clearly as possible. This is not optimized code and the input * matrices have been chosen for simplicity rather than performance. * These should not be used either as a performance guide or for * benchmarking purposes. */ // includes, system #include #include #include #include // CUDA Runtime #include // Using updated (v2) interfaces for CUBLAS and CUSPARSE #include #include // Utilities and system includes #include // shared functions common to CUDA Samples #include // CUDA error checking const char *sSDKname = "conjugateGradientPrecond"; /* * Generate a matrix representing a second order regular Laplacian operator * on a 2D domain in Compressed Sparse Row format. */ void genLaplace(int *row_ptr, int *col_ind, float *val, int M, int N, int nz, float *rhs) { assert(M==N); int n=(int)sqrt((double)N); assert(n*n==N); printf("laplace dimension = %d\n", n); int idx = 0; // loop over degrees of freedom for (int i = 0; i < N; i++) { int ix = i % n; int iy = i / n; row_ptr[i] = idx; // up if (iy > 0) { val[idx] = 1.0; col_ind[idx] = i - n; idx++; } else { rhs[i] -= 1.0; } // left if (ix > 0) { val[idx] = 1.0; col_ind[idx] = i - 1; idx++; } else { rhs[i] -= 0.0; } // center val[idx] = -4.0; col_ind[idx] = i; idx++; //right if (ix < n - 1) { val[idx] = 1.0; col_ind[idx] = i + 1; idx++; } else { rhs[i] -= 0.0; } // down if (iy < n - 1) { val[idx] = 1.0; col_ind[idx] = i + n; idx++; } else { rhs[i] -= 0.0; } } row_ptr[N] = idx; } /* * Solve Ax=b using the conjugate gradient method * a) without any preconditioning, * b) using an Incomplete Cholesky preconditioner, and * c) using an ILU0 preconditioner. */ int main(int argc, char **argv){ const int max_iter = 1000; int k, M = 0, N = 0, nz = 0, *I = NULL, *J = NULL; int *d_col, *d_row; int qatest = 0; const float tol = 1e-12f; float *x, *rhs; float r0, r1, alpha, beta; float *d_val, *d_x; float *d_zm1, *d_zm2, *d_rm2; float *d_r, *d_p, *d_omega, *d_y; float *val = NULL; float *d_valsILU0; float rsum, diff, err = 0.0; float qaerr1, qaerr2 = 0.0; float dot, numerator, denominator, nalpha; const float floatone = 1.0; const float floatzero = 0.0; int nErrors = 0; printf("conjugateGradientPrecond starting...\n"); /* QA testing mode */ if (checkCmdLineFlag(argc, (const char **)argv, "qatest")) { qatest = 1; } /* This will pick the best possible CUDA capable device */ cudaDeviceProp deviceProp; int devID = findCudaDevice(argc, (const char **)argv); printf("GPU selected Device ID = %d \n", devID); if (devID < 0) { printf("Invalid GPU device %d selected, exiting...\n", devID); exit(EXIT_SUCCESS); } checkCudaErrors(cudaGetDeviceProperties(&deviceProp, devID)); /* Statistics about the GPU device */ printf("> GPU device has %d Multi-Processors, SM %d.%d compute capabilities\n\n", deviceProp.multiProcessorCount, deviceProp.major, deviceProp.minor); /* Generate a Laplace matrix in CSR (Compressed Sparse Row) format */ M = N = 16384; nz = 5 * N - 4 * (int)sqrt((double)N); I = (int *)malloc(sizeof(int) * (N + 1)); // csr row pointers for matrix A J = (int *)malloc(sizeof(int) * nz); // csr column indices for matrix A val = (float *)malloc(sizeof(float) * nz); // csr values for matrix A x = (float *)malloc(sizeof(float) * N); rhs = (float *)malloc(sizeof(float) * N); for (int i = 0; i < N; i++) { rhs[i] = 0.0; // Initialize RHS x[i] = 0.0; // Initial solution approximation } genLaplace(I, J, val, M, N, nz, rhs); /* Create CUBLAS context */ cublasHandle_t cublasHandle = NULL; checkCudaErrors(cublasCreate(&cublasHandle)); /* Create CUSPARSE context */ cusparseHandle_t cusparseHandle = NULL; checkCudaErrors(cusparseCreate(&cusparseHandle)); /* Description of the A matrix */ cusparseMatDescr_t descr = 0; checkCudaErrors(cusparseCreateMatDescr(&descr)); checkCudaErrors(cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL)); checkCudaErrors(cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO)); /* Allocate required memory */ checkCudaErrors(cudaMalloc((void **)&d_col, nz * sizeof(int))); checkCudaErrors(cudaMalloc((void **)&d_row, (N + 1) * sizeof(int))); checkCudaErrors(cudaMalloc((void **)&d_val, nz * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_x, N * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_y, N * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_r, N * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_p, N * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_omega, N * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_valsILU0, nz * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_zm1, (N) * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_zm2, (N) * sizeof(float))); checkCudaErrors(cudaMalloc((void **)&d_rm2, (N) * sizeof(float))); /* Wrap raw data into cuSPARSE generic API objects */ cusparseDnVecDescr_t vecp = NULL, vecX=NULL, vecY = NULL, vecR = NULL, vecZM1=NULL; checkCudaErrors(cusparseCreateDnVec(&vecp, N, d_p, CUDA_R_32F)); checkCudaErrors(cusparseCreateDnVec(&vecX, N, d_x, CUDA_R_32F)); checkCudaErrors(cusparseCreateDnVec(&vecY, N, d_y, CUDA_R_32F)); checkCudaErrors(cusparseCreateDnVec(&vecR, N, d_r, CUDA_R_32F)); checkCudaErrors(cusparseCreateDnVec(&vecZM1, N, d_zm1, CUDA_R_32F)); cusparseDnVecDescr_t vecomega = NULL; checkCudaErrors(cusparseCreateDnVec(&vecomega, N, d_omega, CUDA_R_32F)); /* Initialize problem data */ checkCudaErrors(cudaMemcpy( d_col, J, nz * sizeof(int), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_row, I, (N + 1) * sizeof(int), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_val, val, nz * sizeof(float), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_val, val, nz * sizeof(float), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_x, x, N*sizeof(float), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_r, rhs, N * sizeof(float), cudaMemcpyHostToDevice)); cusparseSpMatDescr_t matA = NULL; cusparseSpMatDescr_t matM_lower, matM_upper; cusparseFillMode_t fill_lower = CUSPARSE_FILL_MODE_LOWER; cusparseDiagType_t diag_unit = CUSPARSE_DIAG_TYPE_UNIT; cusparseFillMode_t fill_upper = CUSPARSE_FILL_MODE_UPPER; cusparseDiagType_t diag_non_unit = CUSPARSE_DIAG_TYPE_NON_UNIT; checkCudaErrors(cusparseCreateCsr( &matA, N, N, nz, d_row, d_col, d_val, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_32F)); /* Copy A data to ILU(0) vals as input*/ checkCudaErrors(cudaMemcpy( d_valsILU0, d_val, nz*sizeof(float), cudaMemcpyDeviceToDevice)); //Lower Part checkCudaErrors( cusparseCreateCsr(&matM_lower, N, N, nz, d_row, d_col, d_valsILU0, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_32F) ); checkCudaErrors( cusparseSpMatSetAttribute(matM_lower, CUSPARSE_SPMAT_FILL_MODE, &fill_lower, sizeof(fill_lower)) ); checkCudaErrors( cusparseSpMatSetAttribute(matM_lower, CUSPARSE_SPMAT_DIAG_TYPE, &diag_unit, sizeof(diag_unit)) ); // M_upper checkCudaErrors( cusparseCreateCsr(&matM_upper, N, N, nz, d_row, d_col, d_valsILU0, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_32F) ); checkCudaErrors( cusparseSpMatSetAttribute(matM_upper, CUSPARSE_SPMAT_FILL_MODE, &fill_upper, sizeof(fill_upper)) ); checkCudaErrors( cusparseSpMatSetAttribute(matM_upper, CUSPARSE_SPMAT_DIAG_TYPE, &diag_non_unit, sizeof(diag_non_unit)) ); /* Create ILU(0) info object */ int bufferSizeLU = 0; size_t bufferSizeMV, bufferSizeL, bufferSizeU; void* d_bufferLU, *d_bufferMV, *d_bufferL, *d_bufferU; cusparseSpSVDescr_t spsvDescrL, spsvDescrU; cusparseMatDescr_t matLU; csrilu02Info_t infoILU = NULL; checkCudaErrors(cusparseCreateCsrilu02Info(&infoILU)); checkCudaErrors( cusparseCreateMatDescr(&matLU) ); checkCudaErrors( cusparseSetMatType(matLU, CUSPARSE_MATRIX_TYPE_GENERAL) ); checkCudaErrors( cusparseSetMatIndexBase(matLU, CUSPARSE_INDEX_BASE_ZERO) ); /* Allocate workspace for cuSPARSE */ checkCudaErrors(cusparseSpMV_bufferSize( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matA, vecp, &floatzero, vecomega, CUDA_R_32F, CUSPARSE_SPMV_ALG_DEFAULT, &bufferSizeMV)); checkCudaErrors( cudaMalloc(&d_bufferMV, bufferSizeMV) ); checkCudaErrors(cusparseScsrilu02_bufferSize( cusparseHandle, N, nz, matLU, d_val, d_row, d_col, infoILU, &bufferSizeLU)); checkCudaErrors( cudaMalloc(&d_bufferLU, bufferSizeLU) ); checkCudaErrors( cusparseSpSV_createDescr(&spsvDescrL) ); checkCudaErrors(cusparseSpSV_bufferSize( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_lower, vecR, vecX, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrL, &bufferSizeL)); checkCudaErrors( cudaMalloc(&d_bufferL, bufferSizeL) ); checkCudaErrors( cusparseSpSV_createDescr(&spsvDescrU) ); checkCudaErrors( cusparseSpSV_bufferSize( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_upper, vecR, vecX, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrU, &bufferSizeU)); checkCudaErrors( cudaMalloc(&d_bufferU, bufferSizeU) ); /* Conjugate gradient without preconditioning. ------------------------------------------ Follows the description by Golub & Van Loan, "Matrix Computations 3rd ed.", Section 10.2.6 */ printf("Convergence of CG without preconditioning: \n"); k = 0; r0 = 0; checkCudaErrors(cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1)); while (r1 > tol * tol && k <= max_iter) { k++; if (k == 1) { checkCudaErrors(cublasScopy(cublasHandle, N, d_r, 1, d_p, 1)); } else { beta = r1 / r0; checkCudaErrors(cublasSscal(cublasHandle, N, &beta, d_p, 1)); checkCudaErrors(cublasSaxpy( cublasHandle, N, &floatone, d_r, 1, d_p, 1)); } checkCudaErrors(cusparseSpMV( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matA, vecp, &floatzero, vecomega, CUDA_R_32F, CUSPARSE_SPMV_ALG_DEFAULT, d_bufferMV)); checkCudaErrors(cublasSdot(cublasHandle, N, d_p, 1, d_omega, 1, &dot)); alpha = r1 / dot; checkCudaErrors(cublasSaxpy(cublasHandle, N, &alpha, d_p, 1, d_x, 1)); nalpha = -alpha; checkCudaErrors(cublasSaxpy( cublasHandle, N, &nalpha, d_omega, 1, d_r, 1)); r0 = r1; checkCudaErrors(cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1)); } printf(" iteration = %3d, residual = %e \n", k, sqrt(r1)); checkCudaErrors(cudaMemcpy( x, d_x, N * sizeof(float), cudaMemcpyDeviceToHost)); /* check result */ err = 0.0; for (int i = 0; i < N; i++) { rsum = 0.0; for (int j = I[i]; j < I[i + 1]; j++) { rsum += val[j] * x[J[j]]; } diff = fabs(rsum - rhs[i]); if (diff > err) { err = diff; } } printf(" Convergence Test: %s \n", (k <= max_iter) ? "OK" : "FAIL"); nErrors += (k > max_iter) ? 1 : 0; qaerr1 = err; if (0) { // output result in matlab-style array int n = (int)sqrt((double)N); printf("a = [ "); for (int iy = 0; iy < n; iy++) { for (int ix = 0; ix < n; ix++) { printf(" %f ", x[iy * n + ix]); } if (iy == n - 1) { printf(" ]"); } printf("\n"); } } /* Preconditioned Conjugate Gradient using ILU. -------------------------------------------- Follows the description by Golub & Van Loan, "Matrix Computations 3rd ed.", Algorithm 10.3.1 */ printf("\nConvergence of CG using ILU(0) preconditioning: \n"); /* Perform analysis for ILU(0) */ checkCudaErrors(cusparseScsrilu02_analysis( cusparseHandle, N, nz, descr, d_valsILU0, d_row, d_col, infoILU, CUSPARSE_SOLVE_POLICY_USE_LEVEL, d_bufferLU)); /* generate the ILU(0) factors */ checkCudaErrors(cusparseScsrilu02( cusparseHandle, N, nz, matLU, d_valsILU0, d_row, d_col, infoILU, CUSPARSE_SOLVE_POLICY_USE_LEVEL, d_bufferLU)); /* perform triangular solve analysis */ checkCudaErrors(cusparseSpSV_analysis( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_lower, vecR, vecX, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrL, d_bufferL)); checkCudaErrors(cusparseSpSV_analysis( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_upper, vecR, vecX, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrU, d_bufferU)); /* reset the initial guess of the solution to zero */ for (int i = 0; i < N; i++) { x[i] = 0.0; } checkCudaErrors(cudaMemcpy( d_r, rhs, N * sizeof(float), cudaMemcpyHostToDevice)); checkCudaErrors(cudaMemcpy( d_x, x, N * sizeof(float), cudaMemcpyHostToDevice)); k = 0; checkCudaErrors(cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1)); while (r1 > tol * tol && k <= max_iter) { // preconditioner application: d_zm1 = U^-1 L^-1 d_r checkCudaErrors(cusparseSpSV_solve(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_lower, vecR, vecY, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrL) ); checkCudaErrors(cusparseSpSV_solve(cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matM_upper, vecY, vecZM1, CUDA_R_32F, CUSPARSE_SPSV_ALG_DEFAULT, spsvDescrU)); k++; if (k == 1) { checkCudaErrors(cublasScopy(cublasHandle, N, d_zm1, 1, d_p, 1)); } else { checkCudaErrors(cublasSdot( cublasHandle, N, d_r, 1, d_zm1, 1, &numerator)); checkCudaErrors(cublasSdot( cublasHandle, N, d_rm2, 1, d_zm2, 1, &denominator)); beta = numerator / denominator; checkCudaErrors(cublasSscal(cublasHandle, N, &beta, d_p, 1)); checkCudaErrors(cublasSaxpy( cublasHandle, N, &floatone, d_zm1, 1, d_p, 1)); } checkCudaErrors(cusparseSpMV( cusparseHandle, CUSPARSE_OPERATION_NON_TRANSPOSE, &floatone, matA, vecp, &floatzero, vecomega, CUDA_R_32F, CUSPARSE_SPMV_ALG_DEFAULT, d_bufferMV)); checkCudaErrors(cublasSdot( cublasHandle, N, d_r, 1, d_zm1, 1, &numerator)); checkCudaErrors(cublasSdot( cublasHandle, N, d_p, 1, d_omega, 1, &denominator)); alpha = numerator / denominator; checkCudaErrors(cublasSaxpy(cublasHandle, N, &alpha, d_p, 1, d_x, 1)); checkCudaErrors(cublasScopy(cublasHandle, N, d_r, 1, d_rm2, 1)); checkCudaErrors(cublasScopy(cublasHandle, N, d_zm1, 1, d_zm2, 1)); nalpha = -alpha; checkCudaErrors(cublasSaxpy( cublasHandle, N, &nalpha, d_omega, 1, d_r, 1)); checkCudaErrors(cublasSdot(cublasHandle, N, d_r, 1, d_r, 1, &r1)); } printf(" iteration = %3d, residual = %e \n", k, sqrt(r1)); checkCudaErrors(cudaMemcpy( x, d_x, N * sizeof(float), cudaMemcpyDeviceToHost)); /* check result */ err = 0.0; for (int i = 0; i < N; i++) { rsum = 0.0; for (int j = I[i]; j < I[i + 1]; j++) { rsum += val[j] * x[J[j]]; } diff = fabs(rsum - rhs[i]); if (diff > err) { err = diff; } } printf(" Convergence Test: %s \n", (k <= max_iter) ? "OK" : "FAIL"); nErrors += (k > max_iter) ? 1 : 0; qaerr2 = err; /* Destroy descriptors */ checkCudaErrors(cusparseDestroyCsrilu02Info(infoILU)); checkCudaErrors(cusparseDestroyMatDescr(matLU)); checkCudaErrors(cusparseSpSV_destroyDescr(spsvDescrL)); checkCudaErrors(cusparseSpSV_destroyDescr(spsvDescrU)); checkCudaErrors(cusparseDestroySpMat(matM_lower)); checkCudaErrors(cusparseDestroySpMat(matM_upper)); checkCudaErrors(cusparseDestroySpMat(matA)); checkCudaErrors(cusparseDestroyDnVec(vecp)); checkCudaErrors(cusparseDestroyDnVec(vecomega)); checkCudaErrors(cusparseDestroyDnVec(vecR)); checkCudaErrors(cusparseDestroyDnVec(vecX)); checkCudaErrors(cusparseDestroyDnVec(vecY)); checkCudaErrors(cusparseDestroyDnVec(vecZM1)); /* Destroy contexts */ checkCudaErrors(cusparseDestroy(cusparseHandle)); checkCudaErrors(cublasDestroy(cublasHandle)); /* Free device memory */ free(I); free(J); free(val); free(x); free(rhs); checkCudaErrors(cudaFree(d_bufferMV)); checkCudaErrors(cudaFree(d_bufferLU)); checkCudaErrors(cudaFree(d_bufferL)); checkCudaErrors(cudaFree(d_bufferU)); checkCudaErrors(cudaFree(d_col)); checkCudaErrors(cudaFree(d_row)); checkCudaErrors(cudaFree(d_val)); checkCudaErrors(cudaFree(d_x)); checkCudaErrors(cudaFree(d_y)); checkCudaErrors(cudaFree(d_r)); checkCudaErrors(cudaFree(d_p)); checkCudaErrors(cudaFree(d_omega)); checkCudaErrors(cudaFree(d_valsILU0)); checkCudaErrors(cudaFree(d_zm1)); checkCudaErrors(cudaFree(d_zm2)); checkCudaErrors(cudaFree(d_rm2)); // cudaDeviceReset causes the driver to clean up all state. While // not mandatory in normal operation, it is good practice. It is also // needed to ensure correct operation when the application is being // profiled. Calling cudaDeviceReset causes all profile data to be // flushed before the application exits cudaDeviceReset(); printf("\n"); printf("Test Summary:\n"); printf(" Counted total of %d errors\n", nErrors); printf(" qaerr1 = %f qaerr2 = %f\n\n", fabs(qaerr1), fabs(qaerr2)); exit((nErrors == 0 &&fabs(qaerr1) < 1e-5 && fabs(qaerr2) < 1e-5 ? EXIT_SUCCESS : EXIT_FAILURE)); }