cs_lsolve.c

#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) ;
}


csi cs_lsolve_mrhs (const cs *L, double *x, csi n_rhs)
{
    csi r, 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++)
    {
        for (r = 0; r < n_rhs; r++)
        {
            csi addr = j * n_rhs + r; 
            x [addr] /= Lx [Lp [j]] ;
        }

        for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
        {
            for (r = 0; r < n_rhs; r++)
            {
                csi addr = j * n_rhs + r;
                x [Li [p] * n_rhs + r] -= Lx [p] * x [addr] ;
                //printf("XPU: %f", x [Li [p]]);
            }
        }
    }
    
    return (1) ;
}