LUx = b:
           Ly = b
           Ux = y

def solve(self, b):
        """
        LUx = b:
           Ly = b
           Ux = y
        """
        #b = self.permute_vec(b)
        y = []
        for bri in range(self.N):  # block row index
            for li in range(self.n):  # local row index
                s = 0.0
                for lci in range(li): # local column index
                    s += self.lu[bri][li, lci]*y[bri*self.n+lci]
                for di in range(1, self.ndiag+1):
                    if bri >= di:
                        # di:th sub diagonal (counted as distance from main diag)
                        ci = self.colbyrow[bri][li]
                        s += self.sub[di-1][(bri-di)*self.n+ci]*y[
                            (bri-di)*self.n + ci]
                y.append(b[bri*self.n+self.rowbycol[bri][li]]-s) # Doolittle: L[i, i] == 1
        x = [0]*len(y)
        for bri in range(self.N-1, -1, -1):
            for li in range(self.n - 1, -1, -1):
                s = 0.0
                for ci in range(li+1, self.n):
                    s += self.lu[bri][li, ci]*x[bri*self.n + ci]
                for di in range(1, self.ndiag+1):
                    if bri < self.N-di:
                        ci = self.colbyrow[bri][li]
                        s += self.sup[di-1][bri*self.n+ci]*x[(bri+di)*self.n + ci]
                x[bri*self.n+li] = (y[bri*self.n + li] - s)/self.lu[bri][li, li]
        return x