Finds solution for the rolling case.  Matrix formulation.
166582
def _rolling_lstsq(x, y):
    """"""
    if x.ndim == 2:
        # Treat everything as 3d and avoid AxisError on .swapaxes(1, 2) below
        # This means an original input of:
        #     array([0., 1., 2., 3., 4., 5., 6.])
        # becomes:
        # array([[[0.],
        #         [1.],
        #         [2.],
        #         [3.]],
        #
        #        [[1.],
        #         [2.],
        #         ...
        x = x[:, :, None]
    elif x.ndim <= 1:
        raise np.AxisError("x should have ndmi >= 2")
    return np.squeeze(
        np.matmul(
            np.linalg.inv(np.matmul(x.swapaxes(1, 2), x)),
            np.matmul(x.swapaxes(1, 2), np.atleast_3d(y)),
        )
    )
928518
@Nonnull
  @ReturnsMutableCopy
  public Matrix solve (@Nonnull final Matrix aMatrix)
  {
    if (aMatrix.getRowDimension () != m_nRows)
      throw new IllegalArgumentException ("Matrix row dimensions must agree.");
    if (!isNonsingular ())
      throw new IllegalStateException ("Matrix is singular.");

    // Copy right hand side with pivoting
    final int nCols = aMatrix.getColumnDimension ();
    final Matrix aNewMatrix = aMatrix.getMatrix (m_aPivot, 0, nCols - 1);
    final double [] [] aNewArray = aNewMatrix.internalGetArray ();

    // Solve L*Y = B(piv,:)
    for (int k = 0; k < m_nCols; k++)
    {
      final double [] aNewk = aNewArray[k];
      for (int i = k + 1; i < m_nCols; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    // Solve U*X = Y;
    for (int k = m_nCols - 1; k >= 0; k--)
    {
      final double [] aLUk = m_aLU[k];
      final double [] aNewk = aNewArray[k];
      for (int j = 0; j < nCols; j++)
        aNewk[j] /= aLUk[k];
      for (int i = 0; i < k; i++)
      {
        final double [] aLUi = m_aLU[i];
        final double [] aNewi = aNewArray[i];
        for (int j = 0; j < nCols; j++)
          aNewi[j] -= aNewk[j] * aLUi[k];
      }
    }
    return aNewMatrix;
  }
1067737
def solve(self, lam):
        ''''''
        s = weighted_graphtf_logit(self.nnodes, self.trials, self.successes, lam,
                                 self.Dk.shape[0], self.Dk.shape[1], self.Dk.nnz,
                                 self.Dk.row.astype('int32'), self.Dk.col.astype('int32'), self.Dk.data.astype('double'),
                                 self.maxsteps, self.converge,
                                 self.beta, self.u)
        self.steps.append(s)
        return self.beta
618891
public Matrix solve (Matrix B) throws SingularityException {
      if (B.getRowCount() != m) {
         throw new IllegalArgumentException("Matrix row dimensions must agree.");
      }
      if (!this.isNonsingular()) {
         throw new SingularityException("Matrix is singular.");
      }

      // Copy right hand side with pivoting
      int nx = B.getColumnCount();
      double[][] X = new double[piv.length][nx];
		for (int i=0; i<piv.length; i++) {
			for (int j=0; j<nx; j++) {
				X[i][j]=B.getQuick(piv[i], j);
			}
		}

      // Solve L*Y = B(piv,:)
      for (int k = 0; k < n; k++) {
         for (int i = k+1; i < n; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      // Solve U*X = Y;
      for (int k = n-1; k >= 0; k--) {
         for (int j = 0; j < nx; j++) {
            X[k][j] /= LU[k][k];
         }
         for (int i = 0; i < k; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*LU[i][k];
            }
         }
      }
      return MatrixFactory.createMatrix(X);
   }
927880
public void solve(DenseMatrix B) {
        if (B.nrows() != L.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", L.nrows(), L.ncols(), B.nrows(), B.ncols()));
        }

        int n = B.nrows();
        int nrhs = B.ncols();

        // Solve L*Y = B;
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = 0; i < k; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(k, i));
                }
                B.div(k, j, L.get(k, k));
            }
        }

        // Solve L'*X = Y;
        for (int k = n - 1; k >= 0; k--) {
            for (int j = 0; j < nrhs; j++) {
                for (int i = k + 1; i < n; i++) {
                    B.sub(k, j, B.get(i, j) * L.get(i, k));
                }
                B.div(k, j, L.get(k, k));
            }
        }
    }
1337929
public void solve(DenseMatrix B) {
        if (B.nrows() != L.nrows()) {
            throw new IllegalArgumentException(String.format("Row dimensions do not agree: A is %d x %d, but B is %d x %d", L.nrows(), L.ncols(), B.nrows(), B.ncols()));
        }

        intW info = new intW(0);
        LAPACK.getInstance().dpotrs(NLMatrix.Lower, L.nrows(), B.ncols(), L.data(), L.ld(), B.data(), B.ld(), info);

        if (info.val < 0) {
            logger.error("LAPACK DPOTRS error code: {}", info.val);
            throw new IllegalArgumentException("LAPACK DPOTRS error code: " + info.val);
        }
    }
89704
function fillMatrix(id, t, l, w, h) {
                for (var y = t; y < t + h;) {
                    for (var x = l; x < l + w;) {
                        matrix[y + '-' + x] = id;
                        ++x > maxX && (maxX = x);
                    }
                    ++y > maxY && (maxY = y);
                }
            }
239907
def Rx_matrix(theta):
    """"""
    return np.array([
        [1, 0, 0],
        [0, np.cos(theta), -np.sin(theta)],
        [0, np.sin(theta), np.cos(theta)]
    ])
315240
def cv_compute(self, b, A, B, C, mK, f, m1, m2):
        '''
        

        '''

        A = np.sum([l * a for l, a in zip(self.lam[b], A)
                    if l is not None], axis=0)
        B = np.sum([l * b for l, b in zip(self.lam[b], B)
                    if l is not None], axis=0)
        W = np.linalg.solve(mK + A + C, f)
        if self.transit_model is None:
            model = np.dot(B, W)
        else:
            w_pld = np.concatenate([l * np.dot(self.X(n, m2).T, W)
                                    for n, l in enumerate(self.lam[b])
                                    if l is not None])
            model = np.dot(np.hstack(
                [self.X(n, m1) for n, l in enumerate(self.lam[b])
                 if l is not None]), w_pld)
        model -= np.nanmedian(model)

        return model
926648
public static FloatMatrix solve(FloatMatrix A, FloatMatrix B) {
		A.assertSquare();
		FloatMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.gesv(A.dup(), ipiv, X);
		return X;
	}
1283103
def get_pij_matrix(t, diag, A, A_inv):
    
    return A.dot(np.diag(np.exp(diag * t))).dot(A_inv)
324694
public Matrix solve(Matrix rhs) {
        if (M != N || rhs.M != N || rhs.N != 1) {
            throw new RuntimeException("Illegal matrix dimensions.");
        }

        // create copies of the data
        Matrix A = new Matrix(this);
        Matrix b = new Matrix(rhs);

        // Gaussian elimination with partial pivoting
        for (int i = 0; i < N; i++) {

            // find pivot row and swap
            int max = i;
            for (int j = i + 1; j < N; j++) {
                if (Math.abs(A.data[j][i]) > Math.abs(A.data[max][i])) {
                    max = j;
                }
            }
            A.swap(i, max);
            b.swap(i, max);

            // singular
            if (A.data[i][i] == 0.0) {
                throw new RuntimeException("Matrix is singular.");
            }

            // pivot within b
            for (int j = i + 1; j < N; j++) {
                b.data[j][0] -= b.data[i][0] * A.data[j][i] / A.data[i][i];
            }

            // pivot within A
            for (int j = i + 1; j < N; j++) {
                double m = A.data[j][i] / A.data[i][i];
                for (int k = i + 1; k < N; k++) {
                    A.data[j][k] -= A.data[i][k] * m;
                }
                A.data[j][i] = 0.0;
            }
        }

        // back substitution
        Matrix x = new Matrix(N, 1);
        for (int j = N - 1; j >= 0; j--) {
            double t = 0.0;
            for (int k = j + 1; k < N; k++) {
                t += A.data[j][k] * x.data[k][0];
            }
            x.data[j][0] = (b.data[j][0] - t) / A.data[j][j];
        }
        return x;

    }
199118
def solve_kkt(Q_LU, d, G, A, S_LU, rx, rs, rz, ry):
    
    nineq, nz, neq, nBatch = get_sizes(G, A)

    invQ_rx = rx.btrisolve(*Q_LU)
    if neq > 0:
        h = torch.cat((invQ_rx.unsqueeze(1).bmm(A.transpose(1, 2)).squeeze(1) - ry,
                       invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz), 1)
    else:
        h = invQ_rx.unsqueeze(1).bmm(G.transpose(1, 2)).squeeze(1) + rs / d - rz

    w = -(h.btrisolve(*S_LU))

    g1 = -rx - w[:, neq:].unsqueeze(1).bmm(G).squeeze(1)
    if neq > 0:
        g1 -= w[:, :neq].unsqueeze(1).bmm(A).squeeze(1)
    g2 = -rs - w[:, neq:]

    dx = g1.btrisolve(*Q_LU)
    ds = g2 / d
    dz = w[:, neq:]
    dy = w[:, :neq] if neq > 0 else None

    return dx, ds, dz, dy
641631
public static void solveU(DMatrixSparseCSC U , double []x )
    {
        final int N = U.numCols;

        int idx1 = U.col_idx[N];
        for (int col = N-1; col >= 0; col--) {
            int idx0 = U.col_idx[col];
            double x_j = x[col] /= U.nz_values[idx1-1];

            for (int i = idx0; i < idx1 - 1; i++) {
                int row = U.nz_rows[i];
                x[row] -= U.nz_values[i] * x_j;
            }

            idx1 = idx0;
        }
    }
827776
def aorb(a, b):
    """"""
    return matrix(np.logical_or(a, b).astype('float'), a.size)
787006
public function solve($B)
    {
        if ($B->getRowDimension() == $this->m) {
            if ($this->isFullRank()) {
                // Copy right hand side
                $nx = $B->getColumnDimension();
                $X  = $B->getArrayCopy();
                // Compute Y = transpose(Q)*B
                for ($k = 0; $k < $this->n; $k++) {
                    for ($j = 0; $j < $nx; $j++) {
                        $s = 0.0;
                        for ($i = $k; $i < $this->m; $i++) {
                            $s += $this->QR[$i][$k] * $X[$i][$j];
                        }
                        $s = -$s/$this->QR[$k][$k];
                        for ($i = $k; $i < $this->m; $i++) {
                            $X[$i][$j] += $s * $this->QR[$i][$k];
                        }
                    }
                }
                // Solve R*X = Y;
                for ($k = $this->n-1; $k >= 0; $k--) {
                    for ($j = 0; $j < $nx; $j++) {
                        $X[$k][$j] /= $this->Rdiag[$k];
                    }
                    for ($i = 0; $i < $k; $i++) {
                        for ($j = 0; $j < $nx; $j++) {
                            $X[$i][$j] -= $X[$k][$j]* $this->QR[$i][$k];
                        }
                    }
                }
                $X = new Matrix($X);
                return ($X->getMatrix(0, $this->n-1, 0, $nx));
            } else {
                trigger_error(MatrixRankException, ERROR);
            }
        } else {
            trigger_error(MatrixDimensionException, ERROR);
        }
    }
979000
def MatrixTriangularSolve(a, rhs, lower, adj):
    
    trans = 0 if not adj else 2

    r = np.empty(rhs.shape).astype(a.dtype)
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = sp.linalg.solve_triangular(a[pos] if not adj else np.conj(a[pos]), rhs[pos],
                                            trans=trans, lower=lower)

    return r,