Solve L^T*X = Y

@param X Solution of L*Y=B
@return X
190026
private double[][] solveLtransposed(double[][] X) {
    for(int k = L.length - 1; k >= 0; k--) {
      final double[] Lk = L[k], Xk = X[k];
      timesEquals(Xk, 1. / Lk[k]);
      for(int i = 0; i < k; i++) {
        minusTimesEquals(X[i], Xk, Lk[i]);
      }
    }
    return X;
  }
170317
public static double pow2(int x)
    {
        if(x > Double.MAX_EXPONENT)
            return Double.POSITIVE_INFINITY;
        if(x < Double.MIN_EXPONENT)
            return 0;
        return longBitsToDouble((x+1023L)<<52);
    }
93854
public DoubleMatrix2D solve(DoubleMatrix2D B) {
		if (B.rows() != dim) {
			log.error("wrong dimension of vector b: expected " + dim + ", actual "	+ B.rows());
			throw new RuntimeException("wrong dimension of vector b: expected " + dim	+ ", actual " + B.rows());
		}

	// 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);
		}

		int nOfColumns = B.columns();

		// Solve LY = B (the same as L.Yc = Bc for each column Yc e Bc)
		final double[][] Y = new double[dim][nOfColumns];
		for (int i = 0; i < diagonalLength; i++) {
			double LII = LData[0][i];
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = BI.getQuick(col) / LII;
			}
		}
		for (int i = diagonalLength; i < dim; i++) {
			double[] LI = LData[i - diagonalLength + 1];
			double LII = LI[i];
			double[] sum = new double[nOfColumns];
			for (int j = 0; j < i; j++) {
				double LIJ = LI[j];
				double[] YJ = Y[j];
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LIJ * YJ[col];
				}
			}
			double[] YI = Y[i];
			DoubleMatrix1D BI = B.viewRow(i);
			for (int col = 0; col < nOfColumns; col++) {
				YI[col] = (BI.getQuick(col) - sum[col]) / LII;
			}
		}

		// Solve L[T].X = Y (the same as L[T].Xc = Yc for each column)
		final DoubleMatrix2D X = F2.make(dim, nOfColumns);
		for (int i = dim - 1; i > diagonalLength - 1; i--) {
			double LII = LData[i - diagonalLength + 1][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > i; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}
		for (int i = diagonalLength - 1; i > -1; i--) {
			double LII = LData[0][i];
			double[] sum = new double[nOfColumns];
			for (int j = dim - 1; j > diagonalLength - 1; j--) {
				double[] LJ = LData[j - diagonalLength + 1];
				double LJI = LJ[i];
				DoubleMatrix1D XJ = X.viewRow(j);
				for (int col = 0; col < nOfColumns; col++) {
					sum[col] += LJI * XJ.getQuick(col);
				}
			}
			DoubleMatrix1D XI = X.viewRow(i);
			double[] YI = Y[i];
			for (int col = 0; col < nOfColumns; col++) {
				XI.setQuick(col, (YI[col] - sum[col]) / LII);
			}
		}

		if (this.rescaler != null) {
			// return ALG.mult(this.U, X);
			return ColtUtils.diagonalMatrixMult(this.U, X);
		} else {
			return X;
		}
	}
34314
public static DoubleMatrix asini(DoubleMatrix x) { 
		/*# mapfct('Math.asin') #*/
//RJPP-BEGIN------------------------------------------------------------
	   for (int i = 0; i < x.length; i++)
	      x.put(i, (double) Math.asin(x.get(i)));
	   return x;
//RJPP-END--------------------------------------------------------------
	}
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;
        }
    }
947061
public static void solveL(DMatrixSparseCSC L , double []x )
    {
        final int N = L.numCols;

        int idx0 = L.col_idx[0];
        for (int col = 0; col < N; col++) {

            int idx1 = L.col_idx[col+1];
            double x_j = x[col] /= L.nz_values[idx0];

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

            idx0 = idx1;
        }
    }
848519
public long getIdx(int x,int y){
		long idx=-1;
		if(x>=0&&x<dim1&&y>=0&&y<dim2)
			idx = dim1*x+y;
		return idx;
	}
620343
function(bb, v)
{
    return (bb.l <= v.x && bb.r >= v.x && bb.b <= v.y && bb.t >= v.y);
}
775333
@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");

        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.data[i*2]   = B.data[indexB];
                a.data[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^T)*b = b - u*(gamma*U^T*b)
            for( int n = 0; n < numCols; n++ ) {
                double []u = QR[n];

                double realVV = u[n*2];
                double imagVV = u[n*2+1];

                u[n*2]   = 1;
                u[n*2+1] = 0;

                QrHelperFunctions_ZDRM.rank1UpdateMultR(a, u,0, gammas[n], 0, n, numRows, temp.data);

                u[n*2]   = realVV;
                u[n*2+1] = imagVV;
            }

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

            // save the results
            for( int i = 0; i < numCols; i++ ) {
                int indexB = (i*BnumCols + colB)*2;
                X.data[indexB]   = a.data[i*2];
                X.data[indexB+1] = a.data[i*2+1];
            }
        }
    }
1349005
public double[] solve(double[] b) {
    if(b.length != m) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    if(!this.isNonsingular()) {
      throw new ArithmeticException(ERR_SINGULAR);
    }
    double[] bc = new double[piv.length];
    for(int i = 0; i < piv.length; i++) {
      bc[i] = b[piv[i]];
    }
    solveInplace(bc);
    return n < bc.length ? Arrays.copyOf(bc, n) : bc;
  }
323365
public static DoubleMatrix atani(DoubleMatrix x) { 
		/*# mapfct('Math.atan') #*/
//RJPP-BEGIN------------------------------------------------------------
	   for (int i = 0; i < x.length; i++)
	      x.put(i, (double) Math.atan(x.get(i)));
	   return x;
//RJPP-END--------------------------------------------------------------
	}
135098
@Pure
	public static Point2d L1_L4(double x, double y) {
		final Point2d ntfLambdaPhi = NTFLambert_NTFLambdaPhi(x, y,
				LAMBERT_1_N,
				LAMBERT_1_C,
				LAMBERT_1_XS,
				LAMBERT_1_YS);
		return NTFLambdaPhi_NTFLambert(
				ntfLambdaPhi.getX(), ntfLambdaPhi.getY(),
				LAMBERT_4_N,
				LAMBERT_4_C,
				LAMBERT_4_XS,
				LAMBERT_4_YS);
	}
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;
  }
1148315
public double[][] solve(double[][] B) {
    if(B.length != m) {
      throw new IllegalArgumentException(ERR_MATRIX_DIMENSIONS);
    }
    double[][] sol = solveInplace(getMatrix(B, piv, 0, B[0].length));
    return n < sol.length ? Arrays.copyOf(sol, n) : sol;
  }
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]);
            }
        }
    }
832941
public static double pow(double x, double y) {
        final double lns[] = new double[2];

        if (y == 0.0) {
            return 1.0;
        }

        if (x != x) { // X is NaN
            return x;
        }


        if (x == 0) {
            long bits = Double.doubleToLongBits(x);
            if ((bits & 0x8000000000000000L) != 0) {
                // -zero
                long yi = (long) y;

                if (y < 0 && y == yi && (yi & 1) == 1) {
                    return Double.NEGATIVE_INFINITY;
                }

                if (y > 0 && y == yi && (yi & 1) == 1) {
                    return -0.0;
                }
            }

            if (y < 0) {
                return Double.POSITIVE_INFINITY;
            }
            if (y > 0) {
                return 0.0;
            }

            return Double.NaN;
        }

        if (x == Double.POSITIVE_INFINITY) {
            if (y != y) { // y is NaN
                return y;
            }
            if (y < 0.0) {
                return 0.0;
            } else {
                return Double.POSITIVE_INFINITY;
            }
        }

        if (y == Double.POSITIVE_INFINITY) {
            if (x * x == 1.0) {
                return Double.NaN;
            }

            if (x * x > 1.0) {
                return Double.POSITIVE_INFINITY;
            } else {
                return 0.0;
            }
        }

        if (x == Double.NEGATIVE_INFINITY) {
            if (y != y) { // y is NaN
                return y;
            }

            if (y < 0) {
                long yi = (long) y;
                if (y == yi && (yi & 1) == 1) {
                    return -0.0;
                }

                return 0.0;
            }

            if (y > 0)  {
                long yi = (long) y;
                if (y == yi && (yi & 1) == 1) {
                    return Double.NEGATIVE_INFINITY;
                }

                return Double.POSITIVE_INFINITY;
            }
        }

        if (y == Double.NEGATIVE_INFINITY) {

            if (x * x == 1.0) {
                return Double.NaN;
            }

            if (x * x < 1.0) {
                return Double.POSITIVE_INFINITY;
            } else {
                return 0.0;
            }
        }

        /* Handle special case x<0 */
        if (x < 0) {
            // y is an even integer in this case
            if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
                return pow(-x, y);
            }

            if (y == (long) y) {
                // If y is an integer
                return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
            } else {
                return Double.NaN;
            }
        }

        /* Split y into ya and yb such that y = ya+yb */
        double ya;
        double yb;
        if (y < 8e298 && y > -8e298) {
            double tmp1 = y * HEX_40000000;
            ya = y + tmp1 - tmp1;
            yb = y - ya;
        } else {
            double tmp1 = y * 9.31322574615478515625E-10;
            double tmp2 = tmp1 * 9.31322574615478515625E-10;
            ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
            yb = y - ya;
        }

        /* Compute ln(x) */
        final double lores = log(x, lns);
        if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
            return lores;
        }

        double lna = lns[0];
        double lnb = lns[1];

        /* resplit lns */
        double tmp1 = lna * HEX_40000000;
        double tmp2 = lna + tmp1 - tmp1;
        lnb += lna - tmp2;
        lna = tmp2;

        // y*ln(x) = (aa+ab)
        final double aa = lna * ya;
        final double ab = lna * yb + lnb * ya + lnb * yb;

        lna = aa+ab;
        lnb = -(lna - aa - ab);

        double z = 1.0 / 120.0;
        z = z * lnb + (1.0 / 24.0);
        z = z * lnb + (1.0 / 6.0);
        z = z * lnb + 0.5;
        z = z * lnb + 1.0;
        z = z * lnb;

        final double result = exp(lna, z, null);
        //result = result + result * z;
        return result;
    }
897702
public static FloatMatrix solveSymmetric(FloatMatrix A, FloatMatrix B) {
		A.assertSquare();
		FloatMatrix X = B.dup();
		int[] ipiv = new int[B.rows];
		SimpleBlas.sysv('U', A.dup(), ipiv, X);
		return X;
	}