docs: 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