public static SEXP dgesv(DoubleVector A, DoubleVector B, double tolerance) {
double anorm;
doubleW rcond = new doubleW(0);
if (!Types.isMatrix(A)) {
throw new EvalException("'a' must be a numeric matrix");
}
if (!Types.isMatrix(B)) {
throw new EvalException("'b' must be a numeric matrix");
}
Vector Adims = A.getAttributes().getDim();
Vector Bdims = B.getAttributes().getDim();
int n = Adims.getElementAsInt(0);
if (n == 0) {
throw new EvalException("'a' is 0-diml");
}
int p = Bdims.getElementAsInt(1);
if (p == 0) {
throw new EvalException("no right-hand side in 'b'");
}
if (Adims.getElementAsInt(1) != n) {
throw new EvalException("'a' (" + n + " x " + Adims.getElementAsInt(1) + ") must be square");
}
if (Bdims.getElementAsInt(0) != n) {
throw new EvalException("'b' (" + Bdims.getElementAsInt(0) + " x " + p + ") must be compatible with 'a' (" + n + " x " + n + ")");
}
int ipiv[] = new int[n];
double avals[] = A.toDoubleArray();
LAPACK lapack = LAPACK.getInstance();
intW info = new intW(0);
double[] result = B.toDoubleArray();
lapack.dgesv(n, p, avals, n, ipiv, result, n, info);
if (info.val < 0) {
throw new EvalException("argument -" + info.val + " of Lapack routine 'dgsv' had invalid value");
}
if (info.val > 0) {
throw new EvalException("Lapack routine dgesv: system is exactly singular");
}
anorm = lapack.dlange("1", n, n, A.toDoubleArray(), n, null);
double[] arrWork = new double[4 * n];
lapack.dgecon("1", n, avals, n, anorm, rcond, arrWork, ipiv, info);
if (rcond.val < tolerance) {
throw new EvalException("system is computationally singular: reciprocal condition number = " + rcond.val);
}
return DoubleArrayVector.unsafe(result, B.getAttributes());
}