Newer
Older
#include "cs.h"
/* solve Lx=b where x and b are dense. x=b on input, solution on output. */
csi cs_lsolve (const cs *L, double *x)
{
csi p, j, n, *Lp, *Li ;
double *Lx ;
if (!CS_CSC (L) || !x) return (0) ; /* check inputs */
n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ;
for (j = 0 ; j < n ; j++)
{
x [j] /= Lx [Lp [j]] ;
for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
{
x [Li [p]] -= Lx [p] * x [j] ;
// /* solve Lx=b where x and b are dense. x=b on input, solution on output. */
// //csi cs_lsolve_gpu (const cs *L, double *x)
// // typedef struct cs_sparse /* matrix in compressed-column or triplet form */
// // {
// // csi nzmax ; /* maximum number of entries */
// // csi m ; /* number of rows */
// // csi n ; /* number of columns */
// // csi *p ; /* column pointers (size n+1) or col indices (size nzmax) */
// // csi *i ; /* row indices, size nzmax */
// // double *x ; /* numerical values, size nzmax */
// // csi nz ; /* # of entries in triplet matrix, -1 for compressed-col */
// // } cs ;
// // __global__
// // void cs_lsolve_gpu (int n, const int *Lp, const int *Li, const double *Lx, double *x)
// // {
// // int p, j;
// // // double *Lx ;
// // // n = L->n ;
// // // Lp = L->p ;
// // // Li = L->i ;
// // // Lx = L->x ;
// // for (j = 0 ; j < n ; j++)
// // {
// // x [j] /= Lx [Lp [j]] ;
// // for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
// // {
// // x [Li [p]] -= Lx [p] * x [j] ;
// // }
// // }
// // CUDA kernel to add elements of two arrays
// __global__
// void add(int n, float *x, float *y)
// {
// int index = blockIdx.x * blockDim.x + threadIdx.x;
// int stride = blockDim.x * gridDim.x;
// for (int i = index; i < n; i += stride)
// y[i] = x[i] + y[i];
// }
// // Allocate Unified Memory -- accessible from CPU or GPU
// cudaMallocManaged(&x, N*sizeof(float));
// cudaMallocManaged(&y, N*sizeof(float));
// // initialize x and y arrays on the host
// for (int i = 0; i < N; i++) {
// x[i] = 1.0f;
// y[i] = 2.0f;
// }
// // Prefetch the data to the GPU
// int device = -1;
// cudaGetDevice(&device);
// cudaMemPrefetchAsync(x, N*sizeof(float), device, NULL);
// cudaMemPrefetchAsync(y, N*sizeof(float), device, NULL);
// // Launch kernel on 1M elements on the GPU
// int blockSize = 256;
// int numBlocks = (N + blockSize - 1) / blockSize;
// add<<<numBlocks, blockSize>>>(N, x, y);
// // Wait for GPU to finish before accessing on host
// cudaDeviceSynchronize();
// // Check for errors (all values should be 3.0f)
// float maxError = 0.0f;
// for (int i = 0; i < N; i++)
// maxError = fmax(maxError, fabs(y[i]-3.0f));
// std::cout << "Max error: " << maxError << std::endl;