Skip to content
Snippets Groups Projects
cs_lsolve.c 3.18 KiB
Newer Older
  • Learn to ignore specific revisions
  • Lubomir Riha's avatar
    Lubomir Riha committed
    #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] ;
    
    Lubomir Riha's avatar
    Lubomir Riha committed
                //printf("XPU: %f", x [Li [p]]);
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    
    
    Lubomir Riha's avatar
    Lubomir Riha committed
            }
        }
        return (1) ;
    }
    
    
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // /* solve Lx=b where x and b are dense.  x=b on input, solution on output. */
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // //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;
    
    Lubomir Riha's avatar
    Lubomir Riha committed
        
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // //     // double *Lx ;
    // //     // n = L->n ; 
    // //     // Lp = L->p ; 
    // //     // Li = L->i ; 
    // //     // Lx = L->x ;
    
    Lubomir Riha's avatar
    Lubomir Riha committed
        
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // //     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] ;
    // //         }
    // //     }
    
    Lubomir Riha's avatar
    Lubomir Riha committed
        
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // //     //return (1);
    // // }
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // #include <iostream>
    // #include <math.h>
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // // 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];
    // }
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    // int main(void)
    // {
    //   int N = 1<<20;
    //   float *x, *y;
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // Allocate Unified Memory -- accessible from CPU or GPU
    //   cudaMallocManaged(&x, N*sizeof(float));
    //   cudaMallocManaged(&y, N*sizeof(float));
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // initialize x and y arrays on the host
    //   for (int i = 0; i < N; i++) {
    //     x[i] = 1.0f;
    //     y[i] = 2.0f;
    //   }
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //     // 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);
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // Launch kernel on 1M elements on the GPU
    //   int blockSize = 256;
    //   int numBlocks = (N + blockSize - 1) / blockSize;
    //   add<<<numBlocks, blockSize>>>(N, x, y);
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // Wait for GPU to finish before accessing on host
    //   cudaDeviceSynchronize();
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // 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;
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   // Free memory
    //   cudaFree(x);
    //   cudaFree(y);
    
    Lubomir Riha's avatar
    Lubomir Riha committed
     
    
    Lubomir Riha's avatar
    Lubomir Riha committed
    //   return 0;
    // }