bug: ensure array shape consistency in leastsquares solvers

PyLops solvers output the solution as a nd-shaped array
if the initial model is fed as nd-arrays or if the initial
model is not fed (x0=0) and the operator has forceflat=False. As a side
effect, if x0 is NOT passed as a kwarg (but just as the third arg of
the solver) and the operator forceflat=False, the output is reshaped
even if x0 is a flattened array. I fix this inconsistency making x0
a kwarg when basic solvers are called within a leastsquares solver.

pylops/optimization/cls_leastsquares.py

@@ -248,7 +248,7 @@ def run(
             xinv = cg(
                 self.Op_normal,
                 self.y_normal,
-                self.ncp.zeros(self.Op_normal.shape[1], dtype=self.Op_normal.dtype),
+                x0=self.ncp.zeros(self.Op_normal.shape[1], dtype=self.Op_normal.dtype),
                 **kwargs_solver,
             )[0]
             istop = None
@@ -594,7 +594,7 @@ def run(
             xinv, istop, itn, r1norm, r2norm = cgls(
                 self.RegOp,
                 self.datatot,
-                self.ncp.zeros(self.RegOp.shape[1], dtype=self.RegOp.dtype),
+                x0=self.ncp.zeros(self.RegOp.shape[1], dtype=self.RegOp.dtype),
                 **kwargs_solver,
             )[0:5]
         else:
@@ -816,7 +816,7 @@ def run(
             pinv, istop, itn, r1norm, r2norm = cgls(
                 self.POp,
                 self.y,
-                self.ncp.zeros(self.POp.shape[1], dtype=self.POp.dtype),
+                x0=self.ncp.zeros(self.POp.shape[1], dtype=self.POp.dtype),
                 **kwargs_solver,
             )[0:5]
             # force it 1d as we decorate this method with disable_ndarray_multiplication