#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] ;
            //printf("XPU: %f", x [Li [p]]);

        }
    }
    return (1) ;
}


// /* 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] ;
// //         }
// //     }
    
// //     //return (1);
// // }


    


// #include <iostream>
// #include <math.h>
 
// // 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];
// }
 
// int main(void)
// {
//   int N = 1<<20;
//   float *x, *y;
 
//   // 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;
 
//   // Free memory
//   cudaFree(x);
//   cudaFree(y);
 
//   return 0;
// }