<pre>
-------- CPU linear solver by QR factorization
solve A*x = b, A can be singular
[ls] stands for linear solve
[v] stands for vector
[qr] stands for QR factorization
</pre>
336215
public static int cusolverSpScsrlsvqrHost(
        cusolverSpHandle handle, 
        int m, 
        int nnz, 
        cusparseMatDescr descrA, 
        Pointer csrValA, 
        Pointer csrRowPtrA, 
        Pointer csrColIndA, 
        Pointer b, 
        float tol, 
        int reorder, 
        Pointer x, 
        int[] singularity)
    {
        return checkResult(cusolverSpScsrlsvqrHostNative(handle, m, nnz, descrA, csrValA, csrRowPtrA, csrColIndA, b, tol, reorder, x, singularity));
    }
368902
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X: X rows = "+X.numRows+" expected = "+numCols);
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        U = decomposer.getR(U,true);
        final double gammas[] = decomposer.getGammas();
        final double dataQR[] = QR.data;

        final int BnumCols = B.numCols;

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                a[i*2]   = B.data[indexB];
                a[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            // a = Q_{n-1}...Q_2*Q_1*b
            //
            // Q_n*b = (I-gamma*u*u^H)*b = b - u*(gamma*U^H*b)
            for( int n = 0; n < numCols; n++ ) {
                int indexU = (n*numRows + n + 1)*2;

                double realUb = a[n*2];
                double imagUb = a[n*2+1];

                // U^H*b
                for( int i = n+1; i < numRows; i++ ) {
                    double realU = dataQR[indexU++];
                    double imagU = -dataQR[indexU++];

                    double realB = a[i*2];
                    double imagB = a[i*2+1];

                    realUb += realU*realB - imagU*imagB;
                    imagUb += realU*imagB + imagU*realB;
                }

                // gamma*U^T*b
                realUb *= gammas[n];
                imagUb *= gammas[n];

                a[n*2]   -= realUb;
                a[n*2+1] -= imagUb;

                indexU = (n*numRows + n + 1)*2;
                for( int i = n+1; i < numRows; i++) {
                    double realU = dataQR[indexU++];
                    double imagU = dataQR[indexU++];

                    a[i*2]   -= realU*realUb - imagU*imagUb;
                    a[i*2+1] -= realU*imagUb + imagU*realUb;
                }
            }

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(U.data, a, numCols);

            // save the results

            for( int i = 0; i < numCols; i++ ) {
                int indexX = (i*X.numCols+colB)*2;

                X.data[indexX]   = a[i*2];
                X.data[indexX+1] = a[i*2+1];
            }
        }
    }
785859
private void tql2 (double[][] V) {

   //  This is derived from the Algol procedures tql2, by
   //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
   //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
   //  Fortran subroutine in EISPACK.
   
      for (int i = 1; i < n; i++) {
         e[i-1] = e[i];
      }
      e[n-1] = 0.0;
   
      double f = 0.0;
      double tst1 = 0.0;
      double eps = Math.pow(2.0,-52.0);
      for (int l = 0; l < n; l++) {

         // Find small subdiagonal element
   
         tst1 = Math.max(tst1,Math.abs(d[l]) + Math.abs(e[l]));
         int m = l;
         while (m < n) {
            if (Math.abs(e[m]) <= eps*tst1) {
               break;
            }
            m++;
         }
   
         // If m == l, d[l] is an eigenvalue,
         // otherwise, iterate.
   
         if (m > l) {
            int iter = 0;
            do {
               iter = iter + 1;  // (Could check iteration count here.)
   
               // Compute implicit shift
   
               double g = d[l];
               double p = (d[l+1] - g) / (2.0 * e[l]);
               double r = Math.hypot(p,1.0);
               if (p < 0) {
                  r = -r;
               }
               d[l] = e[l] / (p + r);
               d[l+1] = e[l] * (p + r);
               double dl1 = d[l+1];
               double h = g - d[l];
               for (int i = l+2; i < n; i++) {
                  d[i] -= h;
               }
               f = f + h;
   
               // Implicit QL transformation.
   
               p = d[m];
               double c = 1.0;
               double c2 = c;
               double c3 = c;
               double el1 = e[l+1];
               double s = 0.0;
               double s2 = 0.0;
               for (int i = m-1; i >= l; i--) {
                  c3 = c2;
                  c2 = c;
                  s2 = s;
                  g = c * e[i];
                  h = c * p;
                  r = Math.hypot(p,e[i]);
                  e[i+1] = s * r;
                  s = e[i] / r;
                  c = p / r;
                  p = c * d[i] - s * g;
                  d[i+1] = h + s * (c * g + s * d[i]);
   
                  // Accumulate transformation.
   
                  for (int k = 0; k < n; k++) {
                     h = V[k][i+1];
                     V[k][i+1] = s * V[k][i] + c * h;
                     V[k][i] = c * V[k][i] - s * h;
                  }
               }
               p = -s * s2 * c3 * el1 * e[l] / dl1;
               e[l] = s * p;
               d[l] = c * p;
   
               // Check for convergence.
   
            } while (Math.abs(e[l]) > eps*tst1);
         }
         d[l] = d[l] + f;
         e[l] = 0.0;
      }
     
      // Sort eigenvalues and corresponding vectors.
   
      for (int i = 0; i < n-1; i++) {
         int k = i;
         double p = d[i];
         for (int j = i+1; j < n; j++) {
            if (d[j] < p) {
               k = j;
               p = d[j];
            }
         }
         if (k != i) {
            d[k] = d[i];
            d[i] = p;
            for (int j = 0; j < n; j++) {
               p = V[j][i];
               V[j][i] = V[j][k];
               V[j][k] = p;
            }
         }
      }
   }
223603
@Override
    public void solve(ZMatrixRMaj B, ZMatrixRMaj X) {
        if( X.numRows != numCols )
            throw new IllegalArgumentException("Unexpected dimensions for X");
        else if( B.numRows != numRows || B.numCols != X.numCols )
            throw new IllegalArgumentException("Unexpected dimensions for B");

        int BnumCols = B.numCols;

        Y.reshape(numRows,1);
        Z.reshape(numRows,1);

        // solve each column one by one
        for( int colB = 0; colB < BnumCols; colB++ ) {

            // make a copy of this column in the vector
            for( int i = 0; i < numRows; i++ ) {
                int indexB = B.getIndex(i,colB);
                Y.data[i*2]   = B.data[indexB];
                Y.data[i*2+1] = B.data[indexB+1];
            }

            // Solve Qa=b
            // a = Q'b
            CommonOps_ZDRM.mult(Qt, Y, Z);

            // solve for Rx = b using the standard upper triangular solver
            TriangularSolver_ZDRM.solveU(R.data, Z.data, numCols);

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                X.set(i,colB,Z.data[i*2],Z.data[i*2+1]);
            }
        }
    }
376202
def linreg_ols_qr(y, X):
    """"""
    import numpy as np
    try:  # multiply with inverse to compute coefficients
        q, r = np.linalg.qr(np.dot(X.T, X))
        return np.dot(np.dot(np.linalg.inv(r), q.T), np.dot(X.T, y))
    except np.linalg.LinAlgError:
        print("LinAlgError: Factoring failed")
        return None
500024
public DoubleMatrix1D solve(DoubleMatrix1D b) {
		if (b.size() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ b.size());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + b.size());
		}

		// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
		if (this.rescaler != null) {
			// b = ALG.mult(this.U, b);
			b = ColtUtils.diagonalMatrixMult(this.U, b);
		}

		// Solve L.y = b
		final double[] y = new double[dim];
		for (int i = 0; i < dim; i++) {
			double[] LI = LData[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LI[i];
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve D.z = y
		final double[] z = new double[dim];
		for (int i = 0; i < dim; i++) {
			z[i] = y[i] / DData[i];
		}

		// Solve L[T].x = z
		final DoubleMatrix1D x = F1.make(dim);
		for (int i = dim - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j][i] * x.getQuick(j);
			}
			x.setQuick(i, (z[i] - sum) / LData[i][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}
425543
public void process() {
		solver.solve(Y,X);

		double a = X.data[0];
		double b = X.data[1];
		double c = X.data[2];
		double d = X.data[3];
		double e = X.data[4];

		double bottom = 4*a*c - b*b;
		deltaX = (b*e - 2*d*c)/bottom;
		deltaY = (b*d - 2*a*e)/bottom;
	}
371693
private void _symmetricTql2 ()
  {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.
    for (int i = 1; i < m_nDim; i++)
      m_aEVe[i - 1] = m_aEVe[i];
    m_aEVe[m_nDim - 1] = 0.0;

    double f = 0.0;
    double tst1 = 0.0;
    for (int l = 0; l < m_nDim; l++)
    {
      // Find small subdiagonal element
      tst1 = Math.max (tst1, MathHelper.abs (m_aEVd[l]) + MathHelper.abs (m_aEVe[l]));
      int m = l;
      while (m < m_nDim)
      {
        if (MathHelper.abs (m_aEVe[m]) <= EPSILON * tst1)
          break;
        m++;
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.

      if (m > l)
      {
        do
        {
          // Compute implicit shift
          double g = m_aEVd[l];
          double p = (m_aEVd[l + 1] - g) / (2.0 * m_aEVe[l]);
          double r = MathHelper.hypot (p, 1.0);
          if (p < 0)
          {
            r = -r;
          }
          m_aEVd[l] = m_aEVe[l] / (p + r);
          m_aEVd[l + 1] = m_aEVe[l] * (p + r);
          final double dl1 = m_aEVd[l + 1];
          double h = g - m_aEVd[l];
          for (int i = l + 2; i < m_nDim; i++)
          {
            m_aEVd[i] -= h;
          }
          f += h;

          // Implicit QL transformation.
          p = m_aEVd[m];
          double c = 1.0;
          double c2 = c;
          double c3 = c;
          final double el1 = m_aEVe[l + 1];
          double s = 0.0;
          double s2 = 0.0;
          for (int i = m - 1; i >= l; i--)
          {
            c3 = c2;
            c2 = c;
            s2 = s;
            g = c * m_aEVe[i];
            h = c * p;
            r = MathHelper.hypot (p, m_aEVe[i]);
            m_aEVe[i + 1] = s * r;
            s = m_aEVe[i] / r;
            c = p / r;
            p = c * m_aEVd[i] - s * g;
            m_aEVd[i + 1] = h + s * (c * g + s * m_aEVd[i]);

            // Accumulate transformation.

            for (int k = 0; k < m_nDim; k++)
            {
              h = m_aEigenVector[k][i + 1];
              m_aEigenVector[k][i + 1] = s * m_aEigenVector[k][i] + c * h;
              m_aEigenVector[k][i] = c * m_aEigenVector[k][i] - s * h;
            }
          }
          p = -s * s2 * c3 * el1 * m_aEVe[l] / dl1;
          m_aEVe[l] = s * p;
          m_aEVd[l] = c * p;

          // Check for convergence.

        } while (MathHelper.abs (m_aEVe[l]) > EPSILON * tst1);
      }
      m_aEVd[l] = m_aEVd[l] + f;
      m_aEVe[l] = 0.0;
    }

    // Sort eigenvalues and corresponding vectors.
    for (int i = 0; i < m_nDim - 1; i++)
    {
      int k = i;
      double p = m_aEVd[i];
      for (int j = i + 1; j < m_nDim; j++)
      {
        if (m_aEVd[j] < p)
        {
          k = j;
          p = m_aEVd[j];
        }
      }
      if (k != i)
      {
        m_aEVd[k] = m_aEVd[i];
        m_aEVd[i] = p;
        for (int j = 0; j < m_nDim; j++)
        {
          p = m_aEigenVector[j][i];
          m_aEigenVector[j][i] = m_aEigenVector[j][k];
          m_aEigenVector[j][k] = p;
        }
      }
    }
  }
112739
public static void solveTranL(DMatrixSparseCSC L , double []x )
    {
        final int N = L.numCols;

        for (int j = N-1; j >= 0; j--) {
            int idx0 = L.col_idx[j];
            int idx1 = L.col_idx[j+1];

            for (int p = idx0+1; p < idx1; p++) {
                x[j] -= L.nz_values[p]*x[L.nz_rows[p]];
            }
            x[j] /= L.nz_values[idx0];
        }
    }
41425
def MatrixSolveLs(a, rhs, l2_reg):
    
    r = np.empty(rhs.shape).astype(a.dtype)
    for coord in np.ndindex(a.shape[:-2]):
        pos = coord + (Ellipsis,)
        r[pos] = np.linalg.lstsq(a[pos], rhs[pos])[0]

    return r,
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);
   }
871488
public DoubleMatrix1D solve(DoubleMatrix1D b) {
		if (b.size() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ b.size());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + b.size());
		}

	// with scaling, we must solve U.Q.U.z = U.b, after that we have x = U.z
		if (this.rescaler != null) {
			// b = ALG.mult(this.U, b);
			b = ColtUtils.diagonalMatrixMult(this.U, b);
		}

		// Solve L.y = b
		final double[] y = new double[dim];
		// L is upper left diagonal
		for (int i = 0; i < diagonalLength; i++) {
			double LII = LData[0][i];
			y[i] = b.getQuick(i) / LII;
		}
		for (int i = diagonalLength; i < dim; i++) {
			double[] LI = LData[i - diagonalLength + 1];
			double LII = LI[i];
			double sum = 0;
			for (int j = 0; j < i; j++) {
				sum += LI[j] * y[j];
			}
			y[i] = (b.getQuick(i) - sum) / LII;
		}

		// logger.debug("b: " + ArrayUtils.toString(b));
		// logger.debug("L.y: " + ArrayUtils.toString(getL().operate(y)));

		// Solve L[T].x = y
		final DoubleMatrix1D x = F1.make(dim);
		// for (int i = dim-1; i > -1; i--) {
		// double sum = 0;
		// int ll = Math.max(i, diagonalLength-1);
		// for (int j = dim-1; j > ll; j--) {
		// sum += LData[j][i] * x.getQuick(j);
		//
		// }
		// x.setQuick(i, (y[i] - sum)/ LData[i][i]);
		// }
		for (int i = dim - 1; i > diagonalLength - 1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > i; j--) {
				sum += LData[j - diagonalLength + 1][i] * x.getQuick(j);

			}
			x.setQuick(i, (y[i] - sum) / LData[i - diagonalLength + 1][i]);
		}
		for (int i = diagonalLength - 1; i > -1; i--) {
			double sum = 0;
			for (int j = dim - 1; j > diagonalLength - 1; j--) {
				sum += LData[j - diagonalLength + 1][i] * x.getQuick(j);

			}
			x.setQuick(i, (y[i] - sum) / LData[0][i]);
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, x);
			return ColtUtils.diagonalMatrixMult(this.U, x);
		} else {
			return x;
		}
	}
925836
private double[][] solveInplace(double[][] B) {
    int mx = B.length;
    if(mx != m) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    if(!this.isNonsingular()) {
      throw new ArithmeticException(ERR_SINGULAR);
    }
    // Solve L*Y = B(piv,:)
    for(int k = 0; k < n; k++) {
      final double[] Bk = B[k];
      for(int i = k + 1; i < n; i++) {
        minusTimesEquals(B[i], Bk, LU[i][k]);
      }
    }
    // Solve U*X = Y;
    for(int k = n - 1; k >= 0; k--) {
      final double[] Bk = B[k];
      timesEquals(Bk, 1. / LU[k][k]);
      for(int i = 0; i < k; i++) {
        minusTimesEquals(B[i], Bk, LU[i][k]);
      }
    }
    return B;
  }
903147
def factorize_and_solve(self,A,b):
        

        A = coo_matrix(A)

        x = b.copy()
        self.mumps.set_centralized_assembled_values(A.data)
        self.mumps.set_rhs(x)
        self.mumps.run(job=5)

        return x
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;
  }
792246
protected DMatrixRMaj solveConstraintMatrix() {

		int rowAA = 0;
		for( int i = 0; i < numControl; i++ ) {
			for( int j = i+1; j < numControl; j++ ) {
				for( int k = j; k < numControl; k++ , rowAA++ ) {
					// x_{ii}*x_{jk} = x_{ik}*x_{ji}
					extractXaXb(getIndex(i, i), getIndex(j, k), XiiXjk);
					extractXaXb(getIndex(i, k), getIndex(j, i), XikXji);

					for( int l = 1; l <= AA.numCols; l++ ) {
						AA.set(rowAA,l-1,XikXji[l]-XiiXjk[l]);
					}
					yy.set(rowAA,XiiXjk[0]-XikXji[0]);
				}
			}
		}
//		AA.print();
		CommonOps_DDRM.solve(AA, yy, xx);

		return xx;
	}
921192
@Override
    public ZMatrixRMaj getQ(ZMatrixRMaj Q , boolean compact ) {
        if( compact )
            Q = UtilDecompositons_ZDRM.checkIdentity(Q,numRows,minLength);
        else
            Q = UtilDecompositons_ZDRM.checkIdentity(Q,numRows,numRows);

        for( int j = minLength-1; j >= 0; j-- ) {
            double u[] = dataQR[j];

            double vvReal = u[j*2];
            double vvImag = u[j*2+1];

            u[j*2]   = 1;
            u[j*2+1] = 0;
            double gammaReal = gammas[j];

            QrHelperFunctions_ZDRM.rank1UpdateMultR(Q, u,0, gammaReal,j, j, numRows, v);

            u[j*2]   = vvReal;
            u[j*2+1] = vvImag;
        }

        return Q;
    }