NASP/cuda-samples
85
0
Fork 0
mirror of https://github.com/NVIDIA/cuda-samples.git synced 2026-10-25 22:28:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Apply consistent code formatting across the repo. Add clang-format and pre-commit hooks.

Browse Source
This commit is contained in:
Rob Armstrong
2025-03-27 10:30:07 -07:00
parent 2cd58fbc9a
commit ceab6e8bcc
782 changed files with 107230 additions and 106548 deletions
Download Patch File Download Diff File Expand all files Collapse all files

891
Samples/4_CUDA_Libraries/cuSolverDn_LinearSolver/cuSolverDn_LinearSolver.cpp
View File

@@ -56,529 +56,532 @@
#include <assert.h>
#include <ctype.h>
#include <cuda_runtime.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cuda_runtime.h>
#include "cublas_v2.h"
#include "cusolverDn.h"
#include "helper_cuda.h"
#include "helper_cusolver.h"
template <typename T_ELEM>
int loadMMSparseMatrix(char *filename, char elem_type, bool csrFormat, int *m,
int *n, int *nnz, T_ELEM **aVal, int **aRowInd,
int **aColInd, int extendSymMatrix);
int loadMMSparseMatrix(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
T_ELEM **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
void UsageDN(void) {
printf("<options>\n");
printf("-h : display this help\n");
printf("-R=<name> : choose a linear solver\n");
printf(" chol (cholesky factorization), this is default\n");
printf(" qr (QR factorization)\n");
printf(" lu (LU factorization)\n");
printf("-lda=<int> : leading dimension of A , m by default\n");
printf("-file=<filename>: filename containing a matrix in MM format\n");
printf("-device=<device_id> : <device_id> if want to run on specific GPU\n");
void UsageDN(void)
{
printf("<options>\n");
printf("-h : display this help\n");
printf("-R=<name> : choose a linear solver\n");
printf(" chol (cholesky factorization), this is default\n");
printf(" qr (QR factorization)\n");
printf(" lu (LU factorization)\n");
printf("-lda=<int> : leading dimension of A , m by default\n");
printf("-file=<filename>: filename containing a matrix in MM format\n");
printf("-device=<device_id> : <device_id> if want to run on specific GPU\n");
exit(0);
exit(0);
}
/*
* solve A*x = b by Cholesky factorization
*
*/
int linearSolverCHOL(cusolverDnHandle_t handle, int n, const double *Acopy,
int lda, const double *b, double *x) {
int bufferSize = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
int h_info = 0;
double start, stop;
double time_solve;
cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
int linearSolverCHOL(cusolverDnHandle_t handle, int n, const double *Acopy, int lda, const double *b, double *x)
{
int bufferSize = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
int h_info = 0;
double start, stop;
double time_solve;
cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
checkCudaErrors(cusolverDnDpotrf_bufferSize(handle, uplo, n, (double *)Acopy,
lda, &bufferSize));
checkCudaErrors(cusolverDnDpotrf_bufferSize(handle, uplo, n, (double *)Acopy, lda, &bufferSize));
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
// prepare a copy of A because potrf will overwrite A with L
checkCudaErrors(
cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
// prepare a copy of A because potrf will overwrite A with L
checkCudaErrors(cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
start = second();
start = second();
start = second();
start = second();
checkCudaErrors(
cusolverDnDpotrf(handle, uplo, n, A, lda, buffer, bufferSize, info));
checkCudaErrors(cusolverDnDpotrf(handle, uplo, n, A, lda, buffer, bufferSize, info));
checkCudaErrors(
cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
checkCudaErrors(cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
if (0 != h_info) {
fprintf(stderr, "Error: Cholesky factorization failed\n");
}
if (0 != h_info) {
fprintf(stderr, "Error: Cholesky factorization failed\n");
}
checkCudaErrors(
cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cusolverDnDpotrs(handle, uplo, n, 1, A, lda, x, n, info));
checkCudaErrors(cusolverDnDpotrs(handle, uplo, n, 1, A, lda, x, n, info));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
time_solve = stop - start;
fprintf(stdout, "timing: cholesky = %10.6f sec\n", time_solve);
time_solve = stop - start;
fprintf(stdout, "timing: cholesky = %10.6f sec\n", time_solve);
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
return 0;
return 0;
}
/*
* solve A*x = b by LU with partial pivoting
*
*/
int linearSolverLU(cusolverDnHandle_t handle, int n, const double *Acopy,
int lda, const double *b, double *x) {
int bufferSize = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
int *ipiv = NULL; // pivoting sequence
int h_info = 0;
double start, stop;
double time_solve;
int linearSolverLU(cusolverDnHandle_t handle, int n, const double *Acopy, int lda, const double *b, double *x)
{
int bufferSize = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
int *ipiv = NULL; // pivoting sequence
int h_info = 0;
double start, stop;
double time_solve;
checkCudaErrors(cusolverDnDgetrf_bufferSize(handle, n, n, (double *)Acopy,
lda, &bufferSize));
checkCudaErrors(cusolverDnDgetrf_bufferSize(handle, n, n, (double *)Acopy, lda, &bufferSize));
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
checkCudaErrors(cudaMalloc(&ipiv, sizeof(int) * n));
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
checkCudaErrors(cudaMalloc(&ipiv, sizeof(int) * n));
// prepare a copy of A because getrf will overwrite A with L
checkCudaErrors(
cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
// prepare a copy of A because getrf will overwrite A with L
checkCudaErrors(cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
start = second();
start = second();
start = second();
start = second();
checkCudaErrors(cusolverDnDgetrf(handle, n, n, A, lda, buffer, ipiv, info));
checkCudaErrors(
cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
checkCudaErrors(cusolverDnDgetrf(handle, n, n, A, lda, buffer, ipiv, info));
checkCudaErrors(cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
if (0 != h_info) {
fprintf(stderr, "Error: LU factorization failed\n");
}
if (0 != h_info) {
fprintf(stderr, "Error: LU factorization failed\n");
}
checkCudaErrors(
cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(
cusolverDnDgetrs(handle, CUBLAS_OP_N, n, 1, A, lda, ipiv, x, n, info));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
checkCudaErrors(cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cusolverDnDgetrs(handle, CUBLAS_OP_N, n, 1, A, lda, ipiv, x, n, info));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
time_solve = stop - start;
fprintf(stdout, "timing: LU = %10.6f sec\n", time_solve);
time_solve = stop - start;
fprintf(stdout, "timing: LU = %10.6f sec\n", time_solve);
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
if (ipiv) {
checkCudaErrors(cudaFree(ipiv));
}
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
if (ipiv) {
checkCudaErrors(cudaFree(ipiv));
}
return 0;
return 0;
}
/*
* solve A*x = b by QR
*
*/
int linearSolverQR(cusolverDnHandle_t handle, int n, const double *Acopy,
int lda, const double *b, double *x) {
cublasHandle_t cublasHandle = NULL; // used in residual evaluation
int bufferSize = 0;
int bufferSize_geqrf = 0;
int bufferSize_ormqr = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
double *tau = NULL;
int h_info = 0;
double start, stop;
double time_solve;
const double one = 1.0;
int linearSolverQR(cusolverDnHandle_t handle, int n, const double *Acopy, int lda, const double *b, double *x)
{
cublasHandle_t cublasHandle = NULL; // used in residual evaluation
int bufferSize = 0;
int bufferSize_geqrf = 0;
int bufferSize_ormqr = 0;
int *info = NULL;
double *buffer = NULL;
double *A = NULL;
double *tau = NULL;
int h_info = 0;
double start, stop;
double time_solve;
const double one = 1.0;
checkCudaErrors(cublasCreate(&cublasHandle));
checkCudaErrors(cublasCreate(&cublasHandle));
checkCudaErrors(cusolverDnDgeqrf_bufferSize(handle, n, n, (double *)Acopy,
lda, &bufferSize_geqrf));
checkCudaErrors(cusolverDnDormqr_bufferSize(handle, CUBLAS_SIDE_LEFT,
CUBLAS_OP_T, n, 1, n, A, lda,
NULL, x, n, &bufferSize_ormqr));
checkCudaErrors(cusolverDnDgeqrf_bufferSize(handle, n, n, (double *)Acopy, lda, &bufferSize_geqrf));
checkCudaErrors(cusolverDnDormqr_bufferSize(
handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, n, 1, n, A, lda, NULL, x, n, &bufferSize_ormqr));
printf("buffer_geqrf = %d, buffer_ormqr = %d \n", bufferSize_geqrf,
bufferSize_ormqr);
printf("buffer_geqrf = %d, buffer_ormqr = %d \n", bufferSize_geqrf, bufferSize_ormqr);
bufferSize = (bufferSize_geqrf > bufferSize_ormqr) ? bufferSize_geqrf
: bufferSize_ormqr;
bufferSize = (bufferSize_geqrf > bufferSize_ormqr) ? bufferSize_geqrf : bufferSize_ormqr;
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
checkCudaErrors(cudaMalloc((void **)&tau, sizeof(double) * n));
checkCudaErrors(cudaMalloc(&info, sizeof(int)));
checkCudaErrors(cudaMalloc(&buffer, sizeof(double) * bufferSize));
checkCudaErrors(cudaMalloc(&A, sizeof(double) * lda * n));
checkCudaErrors(cudaMalloc((void **)&tau, sizeof(double) * n));
// prepare a copy of A because getrf will overwrite A with L
checkCudaErrors(
cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
// prepare a copy of A because getrf will overwrite A with L
checkCudaErrors(cudaMemcpy(A, Acopy, sizeof(double) * lda * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
checkCudaErrors(cudaMemset(info, 0, sizeof(int)));
start = second();
start = second();
start = second();
start = second();
// compute QR factorization
checkCudaErrors(
cusolverDnDgeqrf(handle, n, n, A, lda, tau, buffer, bufferSize, info));
// compute QR factorization
checkCudaErrors(cusolverDnDgeqrf(handle, n, n, A, lda, tau, buffer, bufferSize, info));
checkCudaErrors(
cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
checkCudaErrors(cudaMemcpy(&h_info, info, sizeof(int), cudaMemcpyDeviceToHost));
if (0 != h_info) {
fprintf(stderr, "Error: LU factorization failed\n");
}
if (0 != h_info) {
fprintf(stderr, "Error: LU factorization failed\n");
}
checkCudaErrors(
cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
checkCudaErrors(cudaMemcpy(x, b, sizeof(double) * n, cudaMemcpyDeviceToDevice));
// compute Q^T*b
checkCudaErrors(cusolverDnDormqr(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, n, 1,
n, A, lda, tau, x, n, buffer, bufferSize,
info));
// compute Q^T*b
checkCudaErrors(
cusolverDnDormqr(handle, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, n, 1, n, A, lda, tau, x, n, buffer, bufferSize, info));
// x = R \ Q^T*b
checkCudaErrors(cublasDtrsm(cublasHandle, CUBLAS_SIDE_LEFT,
CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_N,
CUBLAS_DIAG_NON_UNIT, n, 1, &one, A, lda, x, n));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
// x = R \ Q^T*b
checkCudaErrors(cublasDtrsm(cublasHandle,
CUBLAS_SIDE_LEFT,
CUBLAS_FILL_MODE_UPPER,
CUBLAS_OP_N,
CUBLAS_DIAG_NON_UNIT,
n,
1,
&one,
A,
lda,
x,
n));
checkCudaErrors(cudaDeviceSynchronize());
stop = second();
time_solve = stop - start;
fprintf(stdout, "timing: QR = %10.6f sec\n", time_solve);
time_solve = stop - start;
fprintf(stdout, "timing: QR = %10.6f sec\n", time_solve);
if (cublasHandle) {
checkCudaErrors(cublasDestroy(cublasHandle));
}
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
if (tau) {
checkCudaErrors(cudaFree(tau));
}
if (cublasHandle) {
checkCudaErrors(cublasDestroy(cublasHandle));
}
if (info) {
checkCudaErrors(cudaFree(info));
}
if (buffer) {
checkCudaErrors(cudaFree(buffer));
}
if (A) {
checkCudaErrors(cudaFree(A));
}
if (tau) {
checkCudaErrors(cudaFree(tau));
}
return 0;
return 0;
}
void parseCommandLineArguments(int argc, char *argv[], struct testOpts &opts) {
memset(&opts, 0, sizeof(opts));
void parseCommandLineArguments(int argc, char *argv[], struct testOpts &opts)
{
memset(&opts, 0, sizeof(opts));
if (checkCmdLineFlag(argc, (const char **)argv, "-h")) {
UsageDN();
}
if (checkCmdLineFlag(argc, (const char **)argv, "R")) {
char *solverType = NULL;
getCmdLineArgumentString(argc, (const char **)argv, "R", &solverType);
if (solverType) {
if ((STRCASECMP(solverType, "chol") != 0) &&
(STRCASECMP(solverType, "lu") != 0) &&
(STRCASECMP(solverType, "qr") != 0)) {
printf("\nIncorrect argument passed to -R option\n");
if (checkCmdLineFlag(argc, (const char **)argv, "-h")) {
UsageDN();
} else {
opts.testFunc = solverType;
}
}
}
if (checkCmdLineFlag(argc, (const char **)argv, "file")) {
char *fileName = 0;
getCmdLineArgumentString(argc, (const char **)argv, "file", &fileName);
if (checkCmdLineFlag(argc, (const char **)argv, "R")) {
char *solverType = NULL;
getCmdLineArgumentString(argc, (const char **)argv, "R", &solverType);
if (fileName) {
opts.sparse_mat_filename = fileName;
} else {
printf("\nIncorrect filename passed to -file \n ");
UsageDN();
}
}
if (checkCmdLineFlag(argc, (const char **)argv, "lda")) {
opts.lda = getCmdLineArgumentInt(argc, (const char **)argv, "lda");
}
}
int main(int argc, char *argv[]) {
struct testOpts opts;
cusolverDnHandle_t handle = NULL;
cublasHandle_t cublasHandle = NULL; // used in residual evaluation
cudaStream_t stream = NULL;
int rowsA = 0; // number of rows of A
int colsA = 0; // number of columns of A
int nnzA = 0; // number of nonzeros of A
int baseA = 0; // base index in CSR format
int lda = 0; // leading dimension in dense matrix
// CSR(A) from I/O
int *h_csrRowPtrA = NULL;
int *h_csrColIndA = NULL;
double *h_csrValA = NULL;
double *h_A = NULL; // dense matrix from CSR(A)
double *h_x = NULL; // a copy of d_x
double *h_b = NULL; // b = ones(m,1)
double *h_r = NULL; // r = b - A*x, a copy of d_r
double *d_A = NULL; // a copy of h_A
double *d_x = NULL; // x = A \ b
double *d_b = NULL; // a copy of h_b
double *d_r = NULL; // r = b - A*x
// the constants are used in residual evaluation, r = b - A*x
const double minus_one = -1.0;
const double one = 1.0;
double x_inf = 0.0;
double r_inf = 0.0;
double A_inf = 0.0;
int errors = 0;
parseCommandLineArguments(argc, argv, opts);
if (NULL == opts.testFunc) {
opts.testFunc = "chol"; // By default running Cholesky as NO solver
// selected with -R option.
}
findCudaDevice(argc, (const char **)argv);
printf("step 1: read matrix market format\n");
if (opts.sparse_mat_filename == NULL) {
opts.sparse_mat_filename = sdkFindFilePath("gr_900_900_crg.mtx", argv[0]);
if (opts.sparse_mat_filename != NULL)
printf("Using default input file [%s]\n", opts.sparse_mat_filename);
else
printf("Could not find gr_900_900_crg.mtx\n");
} else {
printf("Using input file [%s]\n", opts.sparse_mat_filename);
}
if (opts.sparse_mat_filename == NULL) {
fprintf(stderr, "Error: input matrix is not provided\n");
return EXIT_FAILURE;
}
if (loadMMSparseMatrix<double>(opts.sparse_mat_filename, 'd', true, &rowsA,
&colsA, &nnzA, &h_csrValA, &h_csrRowPtrA,
&h_csrColIndA, true)) {
exit(EXIT_FAILURE);
}
baseA = h_csrRowPtrA[0]; // baseA = {0,1}
printf("sparse matrix A is %d x %d with %d nonzeros, base=%d\n", rowsA, colsA,
nnzA, baseA);
if (rowsA != colsA) {
fprintf(stderr, "Error: only support square matrix\n");
exit(EXIT_FAILURE);
}
printf("step 2: convert CSR(A) to dense matrix\n");
lda = opts.lda ? opts.lda : rowsA;
if (lda < rowsA) {
fprintf(stderr, "Error: lda must be greater or equal to dimension of A\n");
exit(EXIT_FAILURE);
}
h_A = (double *)malloc(sizeof(double) * lda * colsA);
h_x = (double *)malloc(sizeof(double) * colsA);
h_b = (double *)malloc(sizeof(double) * rowsA);
h_r = (double *)malloc(sizeof(double) * rowsA);
assert(NULL != h_A);
assert(NULL != h_x);
assert(NULL != h_b);
assert(NULL != h_r);
memset(h_A, 0, sizeof(double) * lda * colsA);
for (int row = 0; row < rowsA; row++) {
const int start = h_csrRowPtrA[row] - baseA;
const int end = h_csrRowPtrA[row + 1] - baseA;
for (int colidx = start; colidx < end; colidx++) {
const int col = h_csrColIndA[colidx] - baseA;
const double Areg = h_csrValA[colidx];
h_A[row + col * lda] = Areg;
}
}
printf("step 3: set right hand side vector (b) to 1\n");
for (int row = 0; row < rowsA; row++) {
h_b[row] = 1.0;
}
// verify if A is symmetric or not.
if (0 == strcmp(opts.testFunc, "chol")) {
int issym = 1;
for (int j = 0; j < colsA; j++) {
for (int i = j; i < rowsA; i++) {
double Aij = h_A[i + j * lda];
double Aji = h_A[j + i * lda];
if (Aij != Aji) {
issym = 0;
break;
if (solverType) {
if ((STRCASECMP(solverType, "chol") != 0) && (STRCASECMP(solverType, "lu") != 0)
&& (STRCASECMP(solverType, "qr") != 0)) {
printf("\nIncorrect argument passed to -R option\n");
UsageDN();
}
else {
opts.testFunc = solverType;
}
}
}
}
if (!issym) {
printf("Error: A has no symmetric pattern, please use LU or QR \n");
exit(EXIT_FAILURE);
if (checkCmdLineFlag(argc, (const char **)argv, "file")) {
char *fileName = 0;
getCmdLineArgumentString(argc, (const char **)argv, "file", &fileName);
if (fileName) {
opts.sparse_mat_filename = fileName;
}
else {
printf("\nIncorrect filename passed to -file \n ");
UsageDN();
}
}
}
checkCudaErrors(cusolverDnCreate(&handle));
checkCudaErrors(cublasCreate(&cublasHandle));
checkCudaErrors(cudaStreamCreate(&stream));
checkCudaErrors(cusolverDnSetStream(handle, stream));
checkCudaErrors(cublasSetStream(cublasHandle, stream));
checkCudaErrors(cudaMalloc((void **)&d_A, sizeof(double) * lda * colsA));
checkCudaErrors(cudaMalloc((void **)&d_x, sizeof(double) * colsA));
checkCudaErrors(cudaMalloc((void **)&d_b, sizeof(double) * rowsA));
checkCudaErrors(cudaMalloc((void **)&d_r, sizeof(double) * rowsA));
printf("step 4: prepare data on device\n");
checkCudaErrors(cudaMemcpy(d_A, h_A, sizeof(double) * lda * colsA,
cudaMemcpyHostToDevice));
checkCudaErrors(
cudaMemcpy(d_b, h_b, sizeof(double) * rowsA, cudaMemcpyHostToDevice));
printf("step 5: solve A*x = b \n");
// d_A and d_b are read-only
if (0 == strcmp(opts.testFunc, "chol")) {
linearSolverCHOL(handle, rowsA, d_A, lda, d_b, d_x);
} else if (0 == strcmp(opts.testFunc, "lu")) {
linearSolverLU(handle, rowsA, d_A, lda, d_b, d_x);
} else if (0 == strcmp(opts.testFunc, "qr")) {
linearSolverQR(handle, rowsA, d_A, lda, d_b, d_x);
} else {
fprintf(stderr, "Error: %s is unknown function\n", opts.testFunc);
exit(EXIT_FAILURE);
}
printf("step 6: evaluate residual\n");
checkCudaErrors(
cudaMemcpy(d_r, d_b, sizeof(double) * rowsA, cudaMemcpyDeviceToDevice));
// r = b - A*x
checkCudaErrors(cublasDgemm_v2(cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, rowsA,
1, colsA, &minus_one, d_A, lda, d_x, rowsA,
&one, d_r, rowsA));
checkCudaErrors(
cudaMemcpy(h_x, d_x, sizeof(double) * colsA, cudaMemcpyDeviceToHost));
checkCudaErrors(
cudaMemcpy(h_r, d_r, sizeof(double) * rowsA, cudaMemcpyDeviceToHost));
x_inf = vec_norminf(colsA, h_x);
r_inf = vec_norminf(rowsA, h_r);
A_inf = mat_norminf(rowsA, colsA, h_A, lda);
printf("|b - A*x| = %E \n", r_inf);
printf("|A| = %E \n", A_inf);
printf("|x| = %E \n", x_inf);
printf("|b - A*x|/(|A|*|x|) = %E \n", r_inf / (A_inf * x_inf));
if (handle) {
checkCudaErrors(cusolverDnDestroy(handle));
}
if (cublasHandle) {
checkCudaErrors(cublasDestroy(cublasHandle));
}
if (stream) {
checkCudaErrors(cudaStreamDestroy(stream));
}
if (h_csrValA) {
free(h_csrValA);
}
if (h_csrRowPtrA) {
free(h_csrRowPtrA);
}
if (h_csrColIndA) {
free(h_csrColIndA);
}
if (h_A) {
free(h_A);
}
if (h_x) {
free(h_x);
}
if (h_b) {
free(h_b);
}
if (h_r) {
free(h_r);
}
if (d_A) {
checkCudaErrors(cudaFree(d_A));
}
if (d_x) {
checkCudaErrors(cudaFree(d_x));
}
if (d_b) {
checkCudaErrors(cudaFree(d_b));
}
if (d_r) {
checkCudaErrors(cudaFree(d_r));
}
return 0;
if (checkCmdLineFlag(argc, (const char **)argv, "lda")) {
opts.lda = getCmdLineArgumentInt(argc, (const char **)argv, "lda");
}
}
int main(int argc, char *argv[])
{
struct testOpts opts;
cusolverDnHandle_t handle = NULL;
cublasHandle_t cublasHandle = NULL; // used in residual evaluation
cudaStream_t stream = NULL;
int rowsA = 0; // number of rows of A
int colsA = 0; // number of columns of A
int nnzA = 0; // number of nonzeros of A
int baseA = 0; // base index in CSR format
int lda = 0; // leading dimension in dense matrix
// CSR(A) from I/O
int *h_csrRowPtrA = NULL;
int *h_csrColIndA = NULL;
double *h_csrValA = NULL;
double *h_A = NULL; // dense matrix from CSR(A)
double *h_x = NULL; // a copy of d_x
double *h_b = NULL; // b = ones(m,1)
double *h_r = NULL; // r = b - A*x, a copy of d_r
double *d_A = NULL; // a copy of h_A
double *d_x = NULL; // x = A \ b
double *d_b = NULL; // a copy of h_b
double *d_r = NULL; // r = b - A*x
// the constants are used in residual evaluation, r = b - A*x
const double minus_one = -1.0;
const double one = 1.0;
double x_inf = 0.0;
double r_inf = 0.0;
double A_inf = 0.0;
int errors = 0;
parseCommandLineArguments(argc, argv, opts);
if (NULL == opts.testFunc) {
opts.testFunc = "chol"; // By default running Cholesky as NO solver
// selected with -R option.
}
findCudaDevice(argc, (const char **)argv);
printf("step 1: read matrix market format\n");
if (opts.sparse_mat_filename == NULL) {
opts.sparse_mat_filename = sdkFindFilePath("gr_900_900_crg.mtx", argv[0]);
if (opts.sparse_mat_filename != NULL)
printf("Using default input file [%s]\n", opts.sparse_mat_filename);
else
printf("Could not find gr_900_900_crg.mtx\n");
}
else {
printf("Using input file [%s]\n", opts.sparse_mat_filename);
}
if (opts.sparse_mat_filename == NULL) {
fprintf(stderr, "Error: input matrix is not provided\n");
return EXIT_FAILURE;
}
if (loadMMSparseMatrix<double>(opts.sparse_mat_filename,
'd',
true,
&rowsA,
&colsA,
&nnzA,
&h_csrValA,
&h_csrRowPtrA,
&h_csrColIndA,
true)) {
exit(EXIT_FAILURE);
}
baseA = h_csrRowPtrA[0]; // baseA = {0,1}
printf("sparse matrix A is %d x %d with %d nonzeros, base=%d\n", rowsA, colsA, nnzA, baseA);
if (rowsA != colsA) {
fprintf(stderr, "Error: only support square matrix\n");
exit(EXIT_FAILURE);
}
printf("step 2: convert CSR(A) to dense matrix\n");
lda = opts.lda ? opts.lda : rowsA;
if (lda < rowsA) {
fprintf(stderr, "Error: lda must be greater or equal to dimension of A\n");
exit(EXIT_FAILURE);
}
h_A = (double *)malloc(sizeof(double) * lda * colsA);
h_x = (double *)malloc(sizeof(double) * colsA);
h_b = (double *)malloc(sizeof(double) * rowsA);
h_r = (double *)malloc(sizeof(double) * rowsA);
assert(NULL != h_A);
assert(NULL != h_x);
assert(NULL != h_b);
assert(NULL != h_r);
memset(h_A, 0, sizeof(double) * lda * colsA);
for (int row = 0; row < rowsA; row++) {
const int start = h_csrRowPtrA[row] - baseA;
const int end = h_csrRowPtrA[row + 1] - baseA;
for (int colidx = start; colidx < end; colidx++) {
const int col = h_csrColIndA[colidx] - baseA;
const double Areg = h_csrValA[colidx];
h_A[row + col * lda] = Areg;
}
}
printf("step 3: set right hand side vector (b) to 1\n");
for (int row = 0; row < rowsA; row++) {
h_b[row] = 1.0;
}
// verify if A is symmetric or not.
if (0 == strcmp(opts.testFunc, "chol")) {
int issym = 1;
for (int j = 0; j < colsA; j++) {
for (int i = j; i < rowsA; i++) {
double Aij = h_A[i + j * lda];
double Aji = h_A[j + i * lda];
if (Aij != Aji) {
issym = 0;
break;
}
}
}
if (!issym) {
printf("Error: A has no symmetric pattern, please use LU or QR \n");
exit(EXIT_FAILURE);
}
}
checkCudaErrors(cusolverDnCreate(&handle));
checkCudaErrors(cublasCreate(&cublasHandle));
checkCudaErrors(cudaStreamCreate(&stream));
checkCudaErrors(cusolverDnSetStream(handle, stream));
checkCudaErrors(cublasSetStream(cublasHandle, stream));
checkCudaErrors(cudaMalloc((void **)&d_A, sizeof(double) * lda * colsA));
checkCudaErrors(cudaMalloc((void **)&d_x, sizeof(double) * colsA));
checkCudaErrors(cudaMalloc((void **)&d_b, sizeof(double) * rowsA));
checkCudaErrors(cudaMalloc((void **)&d_r, sizeof(double) * rowsA));
printf("step 4: prepare data on device\n");
checkCudaErrors(cudaMemcpy(d_A, h_A, sizeof(double) * lda * colsA, cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemcpy(d_b, h_b, sizeof(double) * rowsA, cudaMemcpyHostToDevice));
printf("step 5: solve A*x = b \n");
// d_A and d_b are read-only
if (0 == strcmp(opts.testFunc, "chol")) {
linearSolverCHOL(handle, rowsA, d_A, lda, d_b, d_x);
}
else if (0 == strcmp(opts.testFunc, "lu")) {
linearSolverLU(handle, rowsA, d_A, lda, d_b, d_x);
}
else if (0 == strcmp(opts.testFunc, "qr")) {
linearSolverQR(handle, rowsA, d_A, lda, d_b, d_x);
}
else {
fprintf(stderr, "Error: %s is unknown function\n", opts.testFunc);
exit(EXIT_FAILURE);
}
printf("step 6: evaluate residual\n");
checkCudaErrors(cudaMemcpy(d_r, d_b, sizeof(double) * rowsA, cudaMemcpyDeviceToDevice));
// r = b - A*x
checkCudaErrors(cublasDgemm_v2(
cublasHandle, CUBLAS_OP_N, CUBLAS_OP_N, rowsA, 1, colsA, &minus_one, d_A, lda, d_x, rowsA, &one, d_r, rowsA));
checkCudaErrors(cudaMemcpy(h_x, d_x, sizeof(double) * colsA, cudaMemcpyDeviceToHost));
checkCudaErrors(cudaMemcpy(h_r, d_r, sizeof(double) * rowsA, cudaMemcpyDeviceToHost));
x_inf = vec_norminf(colsA, h_x);
r_inf = vec_norminf(rowsA, h_r);
A_inf = mat_norminf(rowsA, colsA, h_A, lda);
printf("|b - A*x| = %E \n", r_inf);
printf("|A| = %E \n", A_inf);
printf("|x| = %E \n", x_inf);
printf("|b - A*x|/(|A|*|x|) = %E \n", r_inf / (A_inf * x_inf));
if (handle) {
checkCudaErrors(cusolverDnDestroy(handle));
}
if (cublasHandle) {
checkCudaErrors(cublasDestroy(cublasHandle));
}
if (stream) {
checkCudaErrors(cudaStreamDestroy(stream));
}
if (h_csrValA) {
free(h_csrValA);
}
if (h_csrRowPtrA) {
free(h_csrRowPtrA);
}
if (h_csrColIndA) {
free(h_csrColIndA);
}
if (h_A) {
free(h_A);
}
if (h_x) {
free(h_x);
}
if (h_b) {
free(h_b);
}
if (h_r) {
free(h_r);
}
if (d_A) {
checkCudaErrors(cudaFree(d_A));
}
if (d_x) {
checkCudaErrors(cudaFree(d_x));
}
if (d_b) {
checkCudaErrors(cudaFree(d_b));
}
if (d_r) {
checkCudaErrors(cudaFree(d_r));
}
return 0;
}

434
Samples/4_CUDA_Libraries/cuSolverDn_LinearSolver/mmio.c
View File

@@ -1,128 +1,129 @@
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
/* avoid Windows warnings (for example: strcpy, fscanf, etc.) */
#if defined(_WIN32)
#if defined(_WIN32)
#define _CRT_SECURE_NO_WARNINGS
#endif
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
// System includes
#include <ctype.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
// Project includes
#include "mmio.h"
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_)
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_)
{
FILE *f;
FILE *f;
MM_typecode matcode;
int M, N, nz;
int i;
double *val;
int *I, *J;
int M, N, nz;
int i;
double *val;
int *I, *J;
if ((f = fopen(fname, "r")) == NULL)
return -1;
if (mm_read_banner(f, &matcode) != 0)
{
return -1;
if (mm_read_banner(f, &matcode) != 0) {
printf("mm_read_unsymetric: Could not process Matrix Market banner ");
printf(" in file [%s]\n", fname);
return -1;
}
if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
mm_is_sparse(matcode)))
{
if (!(mm_is_real(matcode) && mm_is_matrix(matcode) && mm_is_sparse(matcode))) {
fprintf(stderr, "Sorry, this application does not support ");
fprintf(stderr, "Market Market type: [%s]\n",
mm_typecode_to_str(matcode));
fprintf(stderr, "Market Market type: [%s]\n", mm_typecode_to_str(matcode));
return -1;
}
/* find out size of sparse matrix: M, N, nz .... */
if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
{
if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) {
fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
return -1;
}
*M_ = M;
*N_ = N;
*M_ = M;
*N_ = N;
*nz_ = nz;
/* reseve memory for matrices */
I = (int *) malloc(nz * sizeof(int));
J = (int *) malloc(nz * sizeof(int));
val = (double *) malloc(nz * sizeof(double));
I = (int *)malloc(nz * sizeof(int));
J = (int *)malloc(nz * sizeof(int));
val = (double *)malloc(nz * sizeof(double));
*val_ = val;
*I_ = I;
*J_ = J;
*I_ = I;
*J_ = J;
/* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
/* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
/* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
for (i=0; i<nz; i++)
{
for (i = 0; i < nz; i++) {
if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]) != 3) {
return -1;
}
I[i]--; /* adjust from 1-based to 0-based */
I[i]--; /* adjust from 1-based to 0-based */
J[i]--;
}
fclose(f);
return 0;
}
int mm_is_valid(MM_typecode matcode)
{
if (!mm_is_matrix(matcode)) return 0;
if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) ||
mm_is_skew(matcode))) return 0;
if (!mm_is_matrix(matcode))
return 0;
if (mm_is_dense(matcode) && mm_is_pattern(matcode))
return 0;
if (mm_is_real(matcode) && mm_is_hermitian(matcode))
return 0;
if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || mm_is_skew(matcode)))
return 0;
return 1;
}
int mm_read_banner(FILE *f, MM_typecode *matcode)
{
char line[MM_MAX_LINE_LENGTH];
char banner[MM_MAX_TOKEN_LENGTH];
char mtx[MM_MAX_TOKEN_LENGTH];
char crd[MM_MAX_TOKEN_LENGTH];
char data_type[MM_MAX_TOKEN_LENGTH];
char storage_scheme[MM_MAX_TOKEN_LENGTH];
char line[MM_MAX_LINE_LENGTH];
char banner[MM_MAX_TOKEN_LENGTH];
char mtx[MM_MAX_TOKEN_LENGTH];
char crd[MM_MAX_TOKEN_LENGTH];
char data_type[MM_MAX_TOKEN_LENGTH];
char storage_scheme[MM_MAX_TOKEN_LENGTH];
char *p;
mm_clear_typecode(matcode);
mm_clear_typecode(matcode);
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
return MM_PREMATURE_EOF;
if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type,
storage_scheme) != 5)
if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, storage_scheme) != 5)
return MM_PREMATURE_EOF;
for (p=mtx; *p!='\0'; *p=tolower(*p),p++); /* convert to lower case */
for (p=crd; *p!='\0'; *p=tolower(*p),p++);
for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
for (p = mtx; *p != '\0'; *p = tolower(*p), p++)
; /* convert to lower case */
for (p = crd; *p != '\0'; *p = tolower(*p), p++)
;
for (p = data_type; *p != '\0'; *p = tolower(*p), p++)
;
for (p = storage_scheme; *p != '\0'; *p = tolower(*p), p++)
;
/* check for banner */
if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
@@ -130,7 +131,7 @@ int mm_read_banner(FILE *f, MM_typecode *matcode)
/* first field should be "mtx" */
if (strcmp(mtx, MM_MTX_STR) != 0)
return MM_UNSUPPORTED_TYPE;
return MM_UNSUPPORTED_TYPE;
mm_set_matrix(matcode);
@@ -140,46 +141,39 @@ int mm_read_banner(FILE *f, MM_typecode *matcode)
if (strcmp(crd, MM_SPARSE_STR) == 0)
mm_set_sparse(matcode);
else
if (strcmp(crd, MM_DENSE_STR) == 0)
mm_set_dense(matcode);
else if (strcmp(crd, MM_DENSE_STR) == 0)
mm_set_dense(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* third field */
if (strcmp(data_type, MM_REAL_STR) == 0)
mm_set_real(matcode);
else
if (strcmp(data_type, MM_COMPLEX_STR) == 0)
else if (strcmp(data_type, MM_COMPLEX_STR) == 0)
mm_set_complex(matcode);
else
if (strcmp(data_type, MM_PATTERN_STR) == 0)
else if (strcmp(data_type, MM_PATTERN_STR) == 0)
mm_set_pattern(matcode);
else
if (strcmp(data_type, MM_INT_STR) == 0)
else if (strcmp(data_type, MM_INT_STR) == 0)
mm_set_integer(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* fourth field */
if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
mm_set_general(matcode);
else
if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
else if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
mm_set_symmetric(matcode);
else
if (strcmp(storage_scheme, MM_HERM_STR) == 0)
else if (strcmp(storage_scheme, MM_HERM_STR) == 0)
mm_set_hermitian(matcode);
else
if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
else if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
mm_set_skew(matcode);
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
@@ -188,36 +182,34 @@ int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
{
if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
return MM_COULD_NOT_WRITE_FILE;
else
else
return 0;
}
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz)
{
char line[MM_MAX_LINE_LENGTH];
int num_items_read;
int num_items_read;
/* set return null parameter values, in case we exit with errors */
*M = *N = *nz = 0;
/* now continue scanning until you reach the end-of-comments */
do
{
if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
do {
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
return MM_PREMATURE_EOF;
}while (line[0] == '%');
} while (line[0] == '%');
/* line[] is either blank or has M,N, nz */
if (sscanf(line, "%d %d %d", M, N, nz) == 3)
return 0;
else
do
{
num_items_read = fscanf(f, "%d %d %d", M, N, nz);
if (num_items_read == EOF) return MM_PREMATURE_EOF;
}
while (num_items_read != 3);
do {
num_items_read = fscanf(f, "%d %d %d", M, N, nz);
if (num_items_read == EOF)
return MM_PREMATURE_EOF;
} while (num_items_read != 3);
return 0;
}
@@ -226,28 +218,26 @@ int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
{
char line[MM_MAX_LINE_LENGTH];
int num_items_read;
int num_items_read;
/* set return null parameter values, in case we exit with errors */
*M = *N = 0;
/* now continue scanning until you reach the end-of-comments */
do
{
if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
do {
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
return MM_PREMATURE_EOF;
}while (line[0] == '%');
} while (line[0] == '%');
/* line[] is either blank or has M,N, nz */
if (sscanf(line, "%d %d", M, N) == 2)
return 0;
else /* we have a blank line */
do
{
num_items_read = fscanf(f, "%d %d", M, N);
if (num_items_read == EOF) return MM_PREMATURE_EOF;
}
while (num_items_read != 2);
do {
num_items_read = fscanf(f, "%d %d", M, N);
if (num_items_read == EOF)
return MM_PREMATURE_EOF;
} while (num_items_read != 2);
return 0;
}
@@ -256,75 +246,62 @@ int mm_write_mtx_array_size(FILE *f, int M, int N)
{
if (fprintf(f, "%d %d\n", M, N) != 2)
return MM_COULD_NOT_WRITE_FILE;
else
else
return 0;
}
/*-------------------------------------------------------------------------*/
/******************************************************************/
/* use when I[], J[], and val[]J, and val[] are already allocated */
/******************************************************************/
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode)
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
{
int i;
if (mm_is_complex(matcode))
{
for (i=0; i<nz; i++)
if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
!= 4) return MM_PREMATURE_EOF;
if (mm_is_complex(matcode)) {
for (i = 0; i < nz; i++)
if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2 * i], &val[2 * i + 1]) != 4)
return MM_PREMATURE_EOF;
}
else if (mm_is_real(matcode) || mm_is_integer(matcode))
{
for (i=0; i<nz; i++)
{
if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
!= 3) return MM_PREMATURE_EOF;
else if (mm_is_real(matcode) || mm_is_integer(matcode)) {
for (i = 0; i < nz; i++) {
if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]) != 3)
return MM_PREMATURE_EOF;
}
}
else if (mm_is_pattern(matcode))
{
for (i=0; i<nz; i++)
if (fscanf(f, "%d %d", &I[i], &J[i])
!= 2) return MM_PREMATURE_EOF;
else if (mm_is_pattern(matcode)) {
for (i = 0; i < nz; i++)
if (fscanf(f, "%d %d", &I[i], &J[i]) != 2)
return MM_PREMATURE_EOF;
}
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
double *real, double *imag, MM_typecode matcode)
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *imag, MM_typecode matcode)
{
if (mm_is_complex(matcode))
{
if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
!= 4) return MM_PREMATURE_EOF;
if (mm_is_complex(matcode)) {
if (fscanf(f, "%d %d %lg %lg", I, J, real, imag) != 4)
return MM_PREMATURE_EOF;
}
else if (mm_is_real(matcode) || mm_is_integer(matcode))
{
if (fscanf(f, "%d %d %lg\n", I, J, real)
!= 3) return MM_PREMATURE_EOF;
else if (mm_is_real(matcode) || mm_is_integer(matcode)) {
if (fscanf(f, "%d %d %lg\n", I, J, real) != 3)
return MM_PREMATURE_EOF;
}
else if (mm_is_pattern(matcode))
{
if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
else if (mm_is_pattern(matcode)) {
if (fscanf(f, "%d %d", I, J) != 2)
return MM_PREMATURE_EOF;
}
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
@@ -336,84 +313,78 @@ int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
(nz pairs of real/imaginary values)
************************************************************************/
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
double **val, MM_typecode *matcode)
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, double **val, MM_typecode *matcode)
{
int ret_code;
int ret_code;
FILE *f;
if (strcmp(fname, "stdin") == 0) f=stdin;
else
if ((f = fopen(fname, "r")) == NULL)
if (strcmp(fname, "stdin") == 0)
f = stdin;
else if ((f = fopen(fname, "r")) == NULL)
return MM_COULD_NOT_READ_FILE;
if ((ret_code = mm_read_banner(f, matcode)) != 0)
return ret_code;
if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) &&
mm_is_matrix(*matcode)))
if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && mm_is_matrix(*matcode)))
return MM_UNSUPPORTED_TYPE;
if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
return ret_code;
*I = (int *) malloc(*nz * sizeof(int));
*J = (int *) malloc(*nz * sizeof(int));
*I = (int *)malloc(*nz * sizeof(int));
*J = (int *)malloc(*nz * sizeof(int));
*val = NULL;
if (mm_is_complex(*matcode))
{
*val = (double *) malloc(*nz * 2 * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
if (mm_is_complex(*matcode)) {
*val = (double *)malloc(*nz * 2 * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode);
if (ret_code != 0)
return ret_code;
}
else if (mm_is_real(*matcode) || mm_is_integer(*matcode))
{
*val = (double *) malloc(*nz * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
else if (mm_is_real(*matcode) || mm_is_integer(*matcode)) {
*val = (double *)malloc(*nz * sizeof(double));
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode);
if (ret_code != 0)
return ret_code;
}
else if (mm_is_pattern(*matcode))
{
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val,
*matcode);
if (ret_code != 0) return ret_code;
else if (mm_is_pattern(*matcode)) {
ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, *matcode);
if (ret_code != 0)
return ret_code;
}
if (f != stdin) fclose(f);
if (f != stdin)
fclose(f);
return 0;
}
int mm_write_banner(FILE *f, MM_typecode matcode)
{
char *str = mm_typecode_to_str(matcode);
int ret_code;
int ret_code;
ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
free(str);
if (ret_code !=2 )
if (ret_code != 2)
return MM_COULD_NOT_WRITE_FILE;
else
return 0;
}
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode)
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode)
{
FILE *f;
int i;
int i;
if (strcmp(fname, "stdout") == 0)
if (strcmp(fname, "stdout") == 0)
f = stdout;
else
if ((f = fopen(fname, "w")) == NULL)
else if ((f = fopen(fname, "w")) == NULL)
return MM_COULD_NOT_WRITE_FILE;
/* print banner followed by typecode */
fprintf(f, "%s ", MatrixMarketBanner);
fprintf(f, "%s\n", mm_typecode_to_str(matcode));
@@ -423,54 +394,51 @@ int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
/* print values */
if (mm_is_pattern(matcode))
for (i=0; i<nz; i++)
for (i = 0; i < nz; i++)
fprintf(f, "%d %d\n", I[i], J[i]);
else
if (mm_is_integer(matcode))
for (i=0; i<nz; i++)
else if (mm_is_integer(matcode))
for (i = 0; i < nz; i++)
fprintf(f, "%d %d %d\n", I[i], J[i], (int)val[i]);
else
if (mm_is_real(matcode))
for (i=0; i<nz; i++)
else if (mm_is_real(matcode))
for (i = 0; i < nz; i++)
fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
else
if (mm_is_complex(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i],
val[2*i+1]);
else
{
if (f != stdout) fclose(f);
else if (mm_is_complex(matcode))
for (i = 0; i < nz; i++)
fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2 * i], val[2 * i + 1]);
else {
if (f != stdout)
fclose(f);
return MM_UNSUPPORTED_TYPE;
}
if (f !=stdout) fclose(f);
if (f != stdout)
fclose(f);
return 0;
}
/**
* Create a new copy of a string s. mm_strdup() is a common routine, but
* not part of ANSI C, so it is included here. Used by mm_typecode_to_str().
*
*/
* Create a new copy of a string s. mm_strdup() is a common routine, but
* not part of ANSI C, so it is included here. Used by mm_typecode_to_str().
*
*/
static char *mm_strdup(const char *s)
{
size_t len = strlen(s);
char *s2 = (char *) malloc((len+1)*sizeof(char));
return strcpy(s2, s);
size_t len = strlen(s);
char *s2 = (char *)malloc((len + 1) * sizeof(char));
return strcpy(s2, s);
}
char *mm_typecode_to_str(MM_typecode matcode)
char *mm_typecode_to_str(MM_typecode matcode)
{
char buffer[MM_MAX_LINE_LENGTH];
char buffer[MM_MAX_LINE_LENGTH];
char *types[4];
//char *mm_strdup(const char *);
//int error =0;
// char *mm_strdup(const char *);
// int error =0;
/* check for MTX type */
if (mm_is_matrix(matcode))
if (mm_is_matrix(matcode))
types[0] = MM_MTX_STR;
else
return NULL; // error=1;
@@ -478,8 +446,7 @@ char *mm_typecode_to_str(MM_typecode matcode)
/* check for CRD or ARR matrix */
if (mm_is_sparse(matcode))
types[1] = MM_SPARSE_STR;
else
if (mm_is_dense(matcode))
else if (mm_is_dense(matcode))
types[1] = MM_DENSE_STR;
else
return NULL;
@@ -487,14 +454,11 @@ char *mm_typecode_to_str(MM_typecode matcode)
/* check for element data type */
if (mm_is_real(matcode))
types[2] = MM_REAL_STR;
else
if (mm_is_complex(matcode))
else if (mm_is_complex(matcode))
types[2] = MM_COMPLEX_STR;
else
if (mm_is_pattern(matcode))
else if (mm_is_pattern(matcode))
types[2] = MM_PATTERN_STR;
else
if (mm_is_integer(matcode))
else if (mm_is_integer(matcode))
types[2] = MM_INT_STR;
else
return NULL;
@@ -503,19 +467,15 @@ char *mm_typecode_to_str(MM_typecode matcode)
/* check for symmetry type */
if (mm_is_general(matcode))
types[3] = MM_GENERAL_STR;
else
if (mm_is_symmetric(matcode))
else if (mm_is_symmetric(matcode))
types[3] = MM_SYMM_STR;
else
if (mm_is_hermitian(matcode))
else if (mm_is_hermitian(matcode))
types[3] = MM_HERM_STR;
else
if (mm_is_skew(matcode))
else if (mm_is_skew(matcode))
types[3] = MM_SKEW_STR;
else
return NULL;
sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
sprintf(buffer, "%s %s %s %s", types[0], types[1], types[2], types[3]);
return mm_strdup(buffer);
}

181
Samples/4_CUDA_Libraries/cuSolverDn_LinearSolver/mmio.h
View File

@@ -1,141 +1,138 @@
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
#ifndef MM_IO_H
#define MM_IO_H
#include <stdio.h>
#if defined(__cplusplus)
extern "C" {
extern "C"
{
#endif /* __cplusplus */
#define MM_MAX_LINE_LENGTH 1025
#define MatrixMarketBanner "%%MatrixMarket"
#define MM_MAX_LINE_LENGTH 1025
#define MatrixMarketBanner "%%MatrixMarket"
#define MM_MAX_TOKEN_LENGTH 64
typedef char MM_typecode[4];
typedef char MM_typecode[4];
char *mm_typecode_to_str(MM_typecode matcode);
char *mm_typecode_to_str(MM_typecode matcode);
int mm_read_banner(FILE *f, MM_typecode *matcode);
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
int mm_read_mtx_array_size(FILE *f, int *M, int *N);
int mm_read_banner(FILE *f, MM_typecode *matcode);
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz);
int mm_read_mtx_array_size(FILE *f, int *M, int *N);
int mm_write_banner(FILE *f, MM_typecode matcode);
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
int mm_write_mtx_array_size(FILE *f, int M, int N);
int mm_write_banner(FILE *f, MM_typecode matcode);
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz);
int mm_write_mtx_array_size(FILE *f, int M, int N);
/********************* MM_typecode query fucntions ***************************/
/********************* MM_typecode query fucntions ***************************/
#define mm_is_matrix(typecode) ((typecode)[0]=='M')
#define mm_is_matrix(typecode) ((typecode)[0] == 'M')
#define mm_is_sparse(typecode) ((typecode)[1]=='C')
#define mm_is_coordinate(typecode)((typecode)[1]=='C')
#define mm_is_dense(typecode) ((typecode)[1]=='A')
#define mm_is_array(typecode) ((typecode)[1]=='A')
#define mm_is_sparse(typecode) ((typecode)[1] == 'C')
#define mm_is_coordinate(typecode) ((typecode)[1] == 'C')
#define mm_is_dense(typecode) ((typecode)[1] == 'A')
#define mm_is_array(typecode) ((typecode)[1] == 'A')
#define mm_is_complex(typecode) ((typecode)[2]=='C')
#define mm_is_real(typecode) ((typecode)[2]=='R')
#define mm_is_pattern(typecode) ((typecode)[2]=='P')
#define mm_is_integer(typecode) ((typecode)[2]=='I')
#define mm_is_complex(typecode) ((typecode)[2] == 'C')
#define mm_is_real(typecode) ((typecode)[2] == 'R')
#define mm_is_pattern(typecode) ((typecode)[2] == 'P')
#define mm_is_integer(typecode) ((typecode)[2] == 'I')
#define mm_is_symmetric(typecode)((typecode)[3]=='S')
#define mm_is_general(typecode) ((typecode)[3]=='G')
#define mm_is_skew(typecode) ((typecode)[3]=='K')
#define mm_is_hermitian(typecode)((typecode)[3]=='H')
#define mm_is_symmetric(typecode) ((typecode)[3] == 'S')
#define mm_is_general(typecode) ((typecode)[3] == 'G')
#define mm_is_skew(typecode) ((typecode)[3] == 'K')
#define mm_is_hermitian(typecode) ((typecode)[3] == 'H')
int mm_is_valid(MM_typecode matcode); /* too complex for a macro */
int mm_is_valid(MM_typecode matcode); /* too complex for a macro */
/********************* MM_typecode modify fucntions ***************************/
/********************* MM_typecode modify fucntions ***************************/
#define mm_set_matrix(typecode) ((*typecode)[0]='M')
#define mm_set_coordinate(typecode) ((*typecode)[1]='C')
#define mm_set_array(typecode) ((*typecode)[1]='A')
#define mm_set_dense(typecode) mm_set_array(typecode)
#define mm_set_sparse(typecode) mm_set_coordinate(typecode)
#define mm_set_matrix(typecode) ((*typecode)[0] = 'M')
#define mm_set_coordinate(typecode) ((*typecode)[1] = 'C')
#define mm_set_array(typecode) ((*typecode)[1] = 'A')
#define mm_set_dense(typecode) mm_set_array(typecode)
#define mm_set_sparse(typecode) mm_set_coordinate(typecode)
#define mm_set_complex(typecode)((*typecode)[2]='C')
#define mm_set_real(typecode) ((*typecode)[2]='R')
#define mm_set_pattern(typecode)((*typecode)[2]='P')
#define mm_set_integer(typecode)((*typecode)[2]='I')
#define mm_set_complex(typecode) ((*typecode)[2] = 'C')
#define mm_set_real(typecode) ((*typecode)[2] = 'R')
#define mm_set_pattern(typecode) ((*typecode)[2] = 'P')
#define mm_set_integer(typecode) ((*typecode)[2] = 'I')
#define mm_set_symmetric(typecode)((*typecode)[3]='S')
#define mm_set_general(typecode)((*typecode)[3]='G')
#define mm_set_skew(typecode) ((*typecode)[3]='K')
#define mm_set_hermitian(typecode)((*typecode)[3]='H')
#define mm_set_symmetric(typecode) ((*typecode)[3] = 'S')
#define mm_set_general(typecode) ((*typecode)[3] = 'G')
#define mm_set_skew(typecode) ((*typecode)[3] = 'K')
#define mm_set_hermitian(typecode) ((*typecode)[3] = 'H')
#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
(*typecode)[2]=' ',(*typecode)[3]='G')
#define mm_clear_typecode(typecode) ((*typecode)[0] = (*typecode)[1] = (*typecode)[2] = ' ', (*typecode)[3] = 'G')
#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
/********************* Matrix Market error codes ***************************/
/********************* Matrix Market error codes ***************************/
#define MM_COULD_NOT_READ_FILE 11
#define MM_PREMATURE_EOF 12
#define MM_NOT_MTX 13
#define MM_NO_HEADER 14
#define MM_UNSUPPORTED_TYPE 15
#define MM_LINE_TOO_LONG 16
#define MM_COULD_NOT_WRITE_FILE 17
#define MM_COULD_NOT_READ_FILE 11
#define MM_PREMATURE_EOF 12
#define MM_NOT_MTX 13
#define MM_NO_HEADER 14
#define MM_UNSUPPORTED_TYPE 15
#define MM_LINE_TOO_LONG 16
#define MM_COULD_NOT_WRITE_FILE 17
/******************** Matrix Market internal definitions ********************
/******************** Matrix Market internal definitions ********************
MM_matrix_typecode: 4-character sequence
MM_matrix_typecode: 4-character sequence
ojbect sparse/ data storage
dense type scheme
ojbect sparse/ data storage
dense type scheme
string position: [0] [1] [2] [3]
string position: [0] [1] [2] [3]
Matrix typecode: M(atrix) C(oord) R(eal) G(eneral)
A(array) C(omplex) H(ermitian)
P(attern) S(ymmetric)
I(nteger) K(kew)
Matrix typecode: M(atrix) C(oord) R(eal) G(eneral)
A(array) C(omplex) H(ermitian)
P(attern) S(ymmetric)
I(nteger) K(kew)
***********************************************************************/
***********************************************************************/
#define MM_MTX_STR "matrix"
#define MM_ARRAY_STR "array"
#define MM_DENSE_STR "array"
#define MM_COORDINATE_STR "coordinate"
#define MM_SPARSE_STR "coordinate"
#define MM_COMPLEX_STR "complex"
#define MM_REAL_STR "real"
#define MM_INT_STR "integer"
#define MM_GENERAL_STR "general"
#define MM_SYMM_STR "symmetric"
#define MM_HERM_STR "hermitian"
#define MM_SKEW_STR "skew-symmetric"
#define MM_PATTERN_STR "pattern"
#define MM_MTX_STR "matrix"
#define MM_ARRAY_STR "array"
#define MM_DENSE_STR "array"
#define MM_COORDINATE_STR "coordinate"
#define MM_SPARSE_STR "coordinate"
#define MM_COMPLEX_STR "complex"
#define MM_REAL_STR "real"
#define MM_INT_STR "integer"
#define MM_GENERAL_STR "general"
#define MM_SYMM_STR "symmetric"
#define MM_HERM_STR "hermitian"
#define MM_SKEW_STR "skew-symmetric"
#define MM_PATTERN_STR "pattern"
/* high level routines */
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J,
double **val, MM_typecode *matcode);
/* high level routines */
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, double **val, MM_typecode *matcode);
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
double val[], MM_typecode matcode);
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img,
MM_typecode matcode);
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode);
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[], double val[], MM_typecode matcode);
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J, double *real, double *img, MM_typecode matcode);
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
double **val_, int **I_, int **J_);
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_, double **val_, int **I_, int **J_);
#if defined(__cplusplus)
}
#endif /* __cplusplus */
#endif /* __cplusplus */
#endif

558
Samples/4_CUDA_Libraries/cuSolverDn_LinearSolver/mmio_wrapper.cpp
View File

@@ -1,427 +1,375 @@
#include <cusolverDn.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mmio.h"
#include <cusolverDn.h>
/* avoid Windows warnings (for example: strcpy, fscanf, etc.) */
#if defined(_WIN32)
#if defined(_WIN32)
#define _CRT_SECURE_NO_WARNINGS
#endif
/* various __inline__ __device__ function to initialize a T_ELEM */
template <typename T_ELEM> __inline__ T_ELEM cuGet (int );
template <> __inline__ float cuGet<float >(int x)
{
return float(x);
}
template <typename T_ELEM> __inline__ T_ELEM cuGet(int);
template <> __inline__ float cuGet<float>(int x) { return float(x); }
template <> __inline__ double cuGet<double>(int x)
{
return double(x);
}
template <> __inline__ double cuGet<double>(int x) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(int x)
{
return (make_cuComplex( float(x), 0.0f ));
}
template <> __inline__ cuComplex cuGet<cuComplex>(int x) { return (make_cuComplex(float(x), 0.0f)); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(int x)
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(int x) { return (make_cuDoubleComplex(double(x), 0.0)); }
template <typename T_ELEM> __inline__ T_ELEM cuGet(int, int);
template <> __inline__ float cuGet<float>(int x, int y) { return float(x); }
template <> __inline__ double cuGet<double>(int x, int y) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(int x, int y) { return make_cuComplex(float(x), float(y)); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(int x, int y)
{
return (make_cuDoubleComplex( double(x), 0.0 ));
return (make_cuDoubleComplex(double(x), double(y)));
}
template <typename T_ELEM> __inline__ T_ELEM cuGet (int , int );
template <> __inline__ float cuGet<float >(int x, int y)
{
return float(x);
}
template <typename T_ELEM> __inline__ T_ELEM cuGet(float);
template <> __inline__ float cuGet<float>(float x) { return float(x); }
template <> __inline__ double cuGet<double>(int x, int y)
{
return double(x);
}
template <> __inline__ double cuGet<double>(float x) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(int x, int y)
{
return make_cuComplex( float(x), float(y) );
}
template <> __inline__ cuComplex cuGet<cuComplex>(float x) { return (make_cuComplex(float(x), 0.0f)); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(int x, int y)
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(float x)
{
return (make_cuDoubleComplex( double(x), double(y) ));
return (make_cuDoubleComplex(double(x), 0.0));
}
template <typename T_ELEM> __inline__ T_ELEM cuGet (float );
template <> __inline__ float cuGet<float >(float x)
{
return float(x);
}
template <typename T_ELEM> __inline__ T_ELEM cuGet(float, float);
template <> __inline__ float cuGet<float>(float x, float y) { return float(x); }
template <> __inline__ double cuGet<double>(float x)
{
return double(x);
}
template <> __inline__ double cuGet<double>(float x, float y) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(float x)
{
return (make_cuComplex( float(x), 0.0f ));
}
template <> __inline__ cuComplex cuGet<cuComplex>(float x, float y) { return (make_cuComplex(float(x), float(y))); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(float x)
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(float x, float y)
{
return (make_cuDoubleComplex( double(x), 0.0 ));
return (make_cuDoubleComplex(double(x), double(y)));
}
template <typename T_ELEM> __inline__ T_ELEM cuGet (float, float );
template <> __inline__ float cuGet<float >(float x, float y)
{
return float(x);
}
template <typename T_ELEM> __inline__ T_ELEM cuGet(double);
template <> __inline__ float cuGet<float>(double x) { return float(x); }
template <> __inline__ double cuGet<double>(float x, float y)
{
return double(x);
}
template <> __inline__ double cuGet<double>(double x) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(float x, float y)
{
return (make_cuComplex( float(x), float(y) ));
}
template <> __inline__ cuComplex cuGet<cuComplex>(double x) { return (make_cuComplex(float(x), 0.0f)); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(float x, float y)
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(double x)
{
return (make_cuDoubleComplex( double(x), double(y) ));
return (make_cuDoubleComplex(double(x), 0.0));
}
template <typename T_ELEM> __inline__ T_ELEM cuGet (double );
template <> __inline__ float cuGet<float >(double x)
{
return float(x);
}
template <typename T_ELEM> __inline__ T_ELEM cuGet(double, double);
template <> __inline__ float cuGet<float>(double x, double y) { return float(x); }
template <> __inline__ double cuGet<double>(double x)
{
return double(x);
}
template <> __inline__ double cuGet<double>(double x, double y) { return double(x); }
template <> __inline__ cuComplex cuGet<cuComplex>(double x)
{
return (make_cuComplex( float(x), 0.0f ));
}
template <> __inline__ cuComplex cuGet<cuComplex>(double x, double y) { return (make_cuComplex(float(x), float(y))); }
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(double x)
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(double x, double y)
{
return (make_cuDoubleComplex( double(x), 0.0 ));
return (make_cuDoubleComplex(double(x), double(y)));
}
template <typename T_ELEM> __inline__ T_ELEM cuGet (double, double );
template <> __inline__ float cuGet<float >(double x, double y)
{
return float(x);
}
template <> __inline__ double cuGet<double>(double x, double y)
{
return double(x);
}
template <> __inline__ cuComplex cuGet<cuComplex>(double x, double y)
{
return (make_cuComplex( float(x), float(y) ));
}
template <> __inline__ cuDoubleComplex cuGet<cuDoubleComplex>(double x, double y)
{
return (make_cuDoubleComplex( double(x), double(y) ));
}
static void compress_index(
const int *Ind,
int nnz,
int m,
int *Ptr,
int base)
static void compress_index(const int *Ind, int nnz, int m, int *Ptr, int base)
{
int i;
/* initialize everything to zero */
for(i=0; i<m+1; i++){
Ptr[i]=0;
}
for (i = 0; i < m + 1; i++) {
Ptr[i] = 0;
}
/* count elements in every row */
Ptr[0]=base;
for(i=0; i<nnz; i++){
Ptr[Ind[i]+(1-base)]++;
}
Ptr[0] = base;
for (i = 0; i < nnz; i++) {
Ptr[Ind[i] + (1 - base)]++;
}
/* add all the values */
for(i=0; i<m; i++){
Ptr[i+1]+=Ptr[i];
}
for (i = 0; i < m; i++) {
Ptr[i + 1] += Ptr[i];
}
}
struct cooFormat {
int i ;
int j ;
int p ; // permutation
struct cooFormat
{
int i;
int j;
int p; // permutation
};
int cmp_cooFormat_csr( struct cooFormat *s, struct cooFormat *t)
int cmp_cooFormat_csr(struct cooFormat *s, struct cooFormat *t)
{
if ( s->i < t->i ){
return -1 ;
if (s->i < t->i) {
return -1;
}
else if ( s->i > t->i ){
return 1 ;
else if (s->i > t->i) {
return 1;
}
else{
return s->j - t->j ;
else {
return s->j - t->j;
}
}
int cmp_cooFormat_csc( struct cooFormat *s, struct cooFormat *t)
int cmp_cooFormat_csc(struct cooFormat *s, struct cooFormat *t)
{
if ( s->j < t->j ){
return -1 ;
if (s->j < t->j) {
return -1;
}
else if ( s->j > t->j ){
return 1 ;
else if (s->j > t->j) {
return 1;
}
else{
return s->i - t->i ;
else {
return s->i - t->i;
}
}
typedef int (*FUNPTR) (const void*, const void*) ;
typedef int (*FUNPTR2) ( struct cooFormat *s, struct cooFormat *t) ;
typedef int (*FUNPTR)(const void *, const void *);
typedef int (*FUNPTR2)(struct cooFormat *s, struct cooFormat *t);
static FUNPTR2 fptr_array[2] = {
static FUNPTR2 fptr_array[2] = {
cmp_cooFormat_csr,
cmp_cooFormat_csc,
};
static int verify_pattern(
int m,
int nnz,
int *csrRowPtr,
int *csrColInd)
static int verify_pattern(int m, int nnz, int *csrRowPtr, int *csrColInd)
{
int i, col, start, end, base_index;
int error_found = 0;
if (nnz != (csrRowPtr[m] - csrRowPtr[0])){
fprintf(stderr, "Error (nnz check failed): (csrRowPtr[%d]=%d - csrRowPtr[%d]=%d) != (nnz=%d)\n", 0, csrRowPtr[0], m, csrRowPtr[m], nnz);
if (nnz != (csrRowPtr[m] - csrRowPtr[0])) {
fprintf(stderr,
"Error (nnz check failed): (csrRowPtr[%d]=%d - csrRowPtr[%d]=%d) != (nnz=%d)\n",
0,
csrRowPtr[0],
m,
csrRowPtr[m],
nnz);
error_found = 1;
}
base_index = csrRowPtr[0];
if ((0 != base_index) && (1 != base_index)){
if ((0 != base_index) && (1 != base_index)) {
fprintf(stderr, "Error (base index check failed): base index = %d\n", base_index);
error_found = 1;
}
for (i=0; (!error_found) && (i<m); i++){
start = csrRowPtr[i ] - base_index;
end = csrRowPtr[i+1] - base_index;
if (start > end){
fprintf(stderr, "Error (corrupted row): csrRowPtr[%d] (=%d) > csrRowPtr[%d] (=%d)\n", i, start+base_index, i+1, end+base_index);
for (i = 0; (!error_found) && (i < m); i++) {
start = csrRowPtr[i] - base_index;
end = csrRowPtr[i + 1] - base_index;
if (start > end) {
fprintf(stderr,
"Error (corrupted row): csrRowPtr[%d] (=%d) > csrRowPtr[%d] (=%d)\n",
i,
start + base_index,
i + 1,
end + base_index);
error_found = 1;
}
for (col=start; col<end; col++){
if (csrColInd[col] < base_index){
for (col = start; col < end; col++) {
if (csrColInd[col] < base_index) {
fprintf(stderr, "Error (column vs. base index check failed): csrColInd[%d] < %d\n", col, base_index);
error_found = 1;
}
if ((col < (end-1)) && (csrColInd[col] >= csrColInd[col+1])){
fprintf(stderr, "Error (sorting of the column indecis check failed): (csrColInd[%d]=%d) >= (csrColInd[%d]=%d)\n", col, csrColInd[col], col+1, csrColInd[col+1]);
if ((col < (end - 1)) && (csrColInd[col] >= csrColInd[col + 1])) {
fprintf(
stderr,
"Error (sorting of the column indecis check failed): (csrColInd[%d]=%d) >= (csrColInd[%d]=%d)\n",
col,
csrColInd[col],
col + 1,
csrColInd[col + 1]);
error_found = 1;
}
}
}
return error_found ;
return error_found;
}
template <typename T_ELEM>
int loadMMSparseMatrix(
char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
T_ELEM **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix)
int loadMMSparseMatrix(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
T_ELEM **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix)
{
MM_typecode matcode;
double *tempVal;
int *tempRowInd,*tempColInd;
double *tval;
int *trow,*tcol;
int *csrRowPtr, *cscColPtr;
int i,j,error,base,count;
MM_typecode matcode;
double *tempVal;
int *tempRowInd, *tempColInd;
double *tval;
int *trow, *tcol;
int *csrRowPtr, *cscColPtr;
int i, j, error, base, count;
struct cooFormat *work;
/* read the matrix */
/* read the matrix */
error = mm_read_mtx_crd(filename, m, n, nnz, &trow, &tcol, &tval, &matcode);
if (error) {
fprintf(stderr, "!!!! can not open file: '%s'\n", filename);
return 1;
return 1;
}
/* start error checking */
if (mm_is_complex(matcode) && ((elem_type != 'z') && (elem_type != 'c'))) {
fprintf(stderr, "!!!! complex matrix requires type 'z' or 'c'\n");
return 1;
return 1;
}
if (mm_is_dense(matcode) || mm_is_array(matcode) || mm_is_pattern(matcode) /*|| mm_is_integer(matcode)*/){
if (mm_is_dense(matcode) || mm_is_array(matcode) || mm_is_pattern(matcode) /*|| mm_is_integer(matcode)*/) {
fprintf(stderr, "!!!! dense, array, pattern and integer matrices are not supported\n");
return 1;
return 1;
}
/* if necessary symmetrize the pattern (transform from triangular to full) */
if ((extendSymMatrix) && (mm_is_symmetric(matcode) || mm_is_hermitian(matcode) || mm_is_skew(matcode))){
//count number of non-diagonal elements
count=0;
for(i=0; i<(*nnz); i++){
if (trow[i] != tcol[i]){
if ((extendSymMatrix) && (mm_is_symmetric(matcode) || mm_is_hermitian(matcode) || mm_is_skew(matcode))) {
// count number of non-diagonal elements
count = 0;
for (i = 0; i < (*nnz); i++) {
if (trow[i] != tcol[i]) {
count++;
}
}
//allocate space for the symmetrized matrix
tempRowInd = (int *)malloc((*nnz + count) * sizeof(int));
tempColInd = (int *)malloc((*nnz + count) * sizeof(int));
if (mm_is_real(matcode) || mm_is_integer(matcode)){
// allocate space for the symmetrized matrix
tempRowInd = (int *)malloc((*nnz + count) * sizeof(int));
tempColInd = (int *)malloc((*nnz + count) * sizeof(int));
if (mm_is_real(matcode) || mm_is_integer(matcode)) {
tempVal = (double *)malloc((*nnz + count) * sizeof(double));
}
else{
else {
tempVal = (double *)malloc(2 * (*nnz + count) * sizeof(double));
}
//copy the elements regular and transposed locations
for(j=0, i=0; i<(*nnz); i++){
tempRowInd[j]=trow[i];
tempColInd[j]=tcol[i];
if (mm_is_real(matcode) || mm_is_integer(matcode)){
tempVal[j]=tval[i];
// copy the elements regular and transposed locations
for (j = 0, i = 0; i < (*nnz); i++) {
tempRowInd[j] = trow[i];
tempColInd[j] = tcol[i];
if (mm_is_real(matcode) || mm_is_integer(matcode)) {
tempVal[j] = tval[i];
}
else{
tempVal[2*j] =tval[2*i];
tempVal[2*j+1]=tval[2*i+1];
else {
tempVal[2 * j] = tval[2 * i];
tempVal[2 * j + 1] = tval[2 * i + 1];
}
j++;
if (trow[i] != tcol[i]){
tempRowInd[j]=tcol[i];
tempColInd[j]=trow[i];
if (mm_is_real(matcode) || mm_is_integer(matcode)){
if (mm_is_skew(matcode)){
tempVal[j]=-tval[i];
if (trow[i] != tcol[i]) {
tempRowInd[j] = tcol[i];
tempColInd[j] = trow[i];
if (mm_is_real(matcode) || mm_is_integer(matcode)) {
if (mm_is_skew(matcode)) {
tempVal[j] = -tval[i];
}
else{
tempVal[j]= tval[i];
else {
tempVal[j] = tval[i];
}
}
else{
if(mm_is_hermitian(matcode)){
tempVal[2*j] = tval[2*i];
tempVal[2*j+1]=-tval[2*i+1];
else {
if (mm_is_hermitian(matcode)) {
tempVal[2 * j] = tval[2 * i];
tempVal[2 * j + 1] = -tval[2 * i + 1];
}
else{
tempVal[2*j] = tval[2*i];
tempVal[2*j+1]= tval[2*i+1];
else {
tempVal[2 * j] = tval[2 * i];
tempVal[2 * j + 1] = tval[2 * i + 1];
}
}
j++;
}
}
(*nnz)+=count;
//free temporary storage
(*nnz) += count;
// free temporary storage
free(trow);
free(tcol);
free(tval);
free(tval);
}
else{
tempRowInd=trow;
tempColInd=tcol;
tempVal =tval;
else {
tempRowInd = trow;
tempColInd = tcol;
tempVal = tval;
}
// life time of (trow, tcol, tval) is over.
// please use COO format (tempRowInd, tempColInd, tempVal)
// use qsort to sort COO format
work = (struct cooFormat *)malloc(sizeof(struct cooFormat)*(*nnz));
if (NULL == work){
// use qsort to sort COO format
work = (struct cooFormat *)malloc(sizeof(struct cooFormat) * (*nnz));
if (NULL == work) {
fprintf(stderr, "!!!! allocation error, malloc failed\n");
return 1;
}
for(i=0; i<(*nnz); i++){
for (i = 0; i < (*nnz); i++) {
work[i].i = tempRowInd[i];
work[i].j = tempColInd[i];
work[i].p = i; // permutation is identity
}
if (csrFormat){
/* create row-major ordering of indices (sorted by row and within each row by column) */
qsort(work, *nnz, sizeof(struct cooFormat), (FUNPTR)fptr_array[0] );
}else{
/* create column-major ordering of indices (sorted by column and within each column by row) */
qsort(work, *nnz, sizeof(struct cooFormat), (FUNPTR)fptr_array[1] );
if (csrFormat) {
/* create row-major ordering of indices (sorted by row and within each row by column) */
qsort(work, *nnz, sizeof(struct cooFormat), (FUNPTR)fptr_array[0]);
}
else {
/* create column-major ordering of indices (sorted by column and within each column by row) */
qsort(work, *nnz, sizeof(struct cooFormat), (FUNPTR)fptr_array[1]);
}
// (tempRowInd, tempColInd) is sorted either by row-major or by col-major
for(i=0; i<(*nnz); i++){
for (i = 0; i < (*nnz); i++) {
tempRowInd[i] = work[i].i;
tempColInd[i] = work[i].j;
}
// setup base
// setup base
// check if there is any row/col 0, if so base-0
// check if there is any row/col equal to matrix dimension m/n, if so base-1
int base0 = 0;
int base1 = 0;
for(i=0; i<(*nnz); i++){
for (i = 0; i < (*nnz); i++) {
const int row = tempRowInd[i];
const int col = tempColInd[i];
if ( (0 == row) || (0 == col) ){
if ((0 == row) || (0 == col)) {
base0 = 1;
}
if ( (*m == row) || (*n == col) ){
if ((*m == row) || (*n == col)) {
base1 = 1;
}
}
if ( base0 && base1 ){
if (base0 && base1) {
printf("Error: input matrix is base-0 and base-1 \n");
return 1;
}
base = 0;
if (base1){
if (base1) {
base = 1;
}
/* compress the appropriate indices */
if (csrFormat){
if (csrFormat) {
/* CSR format (assuming row-major format) */
csrRowPtr = (int *)malloc(((*m)+1) * sizeof(csrRowPtr[0]));
if (!csrRowPtr) return 1;
csrRowPtr = (int *)malloc(((*m) + 1) * sizeof(csrRowPtr[0]));
if (!csrRowPtr)
return 1;
compress_index(tempRowInd, *nnz, *m, csrRowPtr, base);
*aRowInd = csrRowPtr;
@@ -429,101 +377,97 @@ int loadMMSparseMatrix(
}
else {
/* CSC format (assuming column-major format) */
cscColPtr = (int *)malloc(((*n)+1) * sizeof(cscColPtr[0]));
if (!cscColPtr) return 1;
cscColPtr = (int *)malloc(((*n) + 1) * sizeof(cscColPtr[0]));
if (!cscColPtr)
return 1;
compress_index(tempColInd, *nnz, *n, cscColPtr, base);
*aColInd = cscColPtr;
*aRowInd = (int *)malloc((*nnz) * sizeof(int));
}
}
/* transfrom the matrix values of type double into one of the cusparse library types */
/* transfrom the matrix values of type double into one of the cusparse library types */
*aVal = (T_ELEM *)malloc((*nnz) * sizeof(T_ELEM));
for (i=0; i<(*nnz); i++) {
if (csrFormat){
for (i = 0; i < (*nnz); i++) {
if (csrFormat) {
(*aColInd)[i] = tempColInd[i];
}
else{
else {
(*aRowInd)[i] = tempRowInd[i];
}
if (mm_is_real(matcode) || mm_is_integer(matcode)){
(*aVal)[i] = cuGet<T_ELEM>( tempVal[ work[i].p ] );
if (mm_is_real(matcode) || mm_is_integer(matcode)) {
(*aVal)[i] = cuGet<T_ELEM>(tempVal[work[i].p]);
}
else{
(*aVal)[i] = cuGet<T_ELEM>(tempVal[2*work[i].p], tempVal[2*work[i].p+1]);
else {
(*aVal)[i] = cuGet<T_ELEM>(tempVal[2 * work[i].p], tempVal[2 * work[i].p + 1]);
}
}
/* check for corruption */
int error_found;
if (csrFormat){
if (csrFormat) {
error_found = verify_pattern(*m, *nnz, *aRowInd, *aColInd);
}else{
}
else {
error_found = verify_pattern(*n, *nnz, *aColInd, *aRowInd);
}
if (error_found){
if (error_found) {
fprintf(stderr, "!!!! verify_pattern failed\n");
return 1;
}
/* cleanup and exit */
free(work);
free(tempVal);
free(tempVal);
free(tempColInd);
free(tempRowInd);
return 0;
}
}
/* specific instantiation */
template int loadMMSparseMatrix<float>(
char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
float **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<float>(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
float **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<double>(
char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
double **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<cuComplex>(
char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
cuComplex **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<cuDoubleComplex>(
char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
cuDoubleComplex **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<double>(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
double **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<cuComplex>(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
cuComplex **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);
template int loadMMSparseMatrix<cuDoubleComplex>(char *filename,
char elem_type,
bool csrFormat,
int *m,
int *n,
int *nnz,
cuDoubleComplex **aVal,
int **aRowInd,
int **aColInd,
int extendSymMatrix);