research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Another delta on the dogleg implementation of the nonlinear solver.

Browse source 
Added lapack routines for calculation of condition number and to 
fill out the QR factorization capability.
This commit is contained in:
Harry Moffat 2011-02-19 00:08:29 +00:00
parent 08716fa3f1
commit 144c6b23bd
14 changed files with 2200 additions and 35 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
Cantera/src/base/Array.h
View file

@ -120,6 +120,16 @@ namespace Cantera {
m_data.resize(n*m, v);
}
//! Copy the data from one array into another without doing any checking
/*!
* This differs from the assignment operator as no resizing is done and memcpy() is used.
* @param y Array to be copied
*/
void copyData(const Array2D& y) {
size_t n = sizeof(doublereal) * m_nrows * m_ncols;
(void) memcpy(DATA_PTR(m_data), y.ptrColumn(0), n);
}
//! Append a column to the existing matrix using a std vector
/*!
* This operation will add a column onto the existing matrix.

4
Cantera/src/equil/vcs_solve.h
View file

@ -7,7 +7,7 @@
* $Id$
*/
/*
* Copywrite (2005) Sandia Corporation. Under the terms of
*_ Copywrite (2005) Sandia Corporation. Under the terms of
* Contract DE-AC04-94AL85000 with Sandia Corporation, the
* U.S. Government retains certain rights in this software.
*/
@ -2073,6 +2073,8 @@ public:
*/
int m_VCS_UnitsFormat;
friend class vcs_phaseStabilitySolve;
};
#ifdef ALTLINPROG

173
Cantera/src/numerics/NonlinearSolver.cpp
View file

@ -127,15 +127,22 @@ namespace Cantera {
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopy_(0),
deltax_cp_(0),
deltaX_CP_(0),
deltaX_Newton_(0),
residNorm2Cauchy_(0.0),
RJd_norm_(0.0),
lambda_(0.0),
Jd_(0),
trustDeltaX_(0),
trustDelta_(1.0)
deltaX_trust_(0),
trustDelta_(1.0),
Nuu_(0.0),
dist_R0_(0.0),
dist_R1_(0.0),
dist_R2_(0.0),
dist_Total_(0.0),
JdJd_norm_(0.0),
normTrust_Newton_(0.0),
normTrust_CP_(0.0)
{
neq_ = m_func->nEquations();
@ -164,9 +171,9 @@ namespace Cantera {
#ifdef DEBUG_DOGLEG
jacCopy_.resize(neq_, neq_, 0.0);
deltax_cp_.resize(neq_, 0.0);
deltaX_CP_.resize(neq_, 0.0);
Jd_.resize(neq_, 0.0);
trustDeltaX_.resize(neq_, 1.0);
deltaX_trust_.resize(neq_, 1.0);
#endif
}
@ -217,14 +224,22 @@ namespace Cantera {
m_print_flag(0),
m_ScaleSolnNormToResNorm(0.001),
jacCopy_(0),
deltax_cp_(0),
deltaX_CP_(0),
deltaX_Newton_(0),
residNorm2Cauchy_(0.0),
RJd_norm_(0.0),
lambda_(0.0),
Jd_(0),
trustDeltaX_(0),
trustDelta_(1.0)
deltaX_trust_(0),
trustDelta_(1.0),
Nuu_(0.0),
dist_R0_(0.0),
dist_R1_(0.0),
dist_R2_(0.0),
dist_Total_(0.0),
JdJd_norm_(0.0),
normTrust_Newton_(0.0),
normTrust_CP_(0.0)
{
*this =operator=(right);
}
@ -285,14 +300,23 @@ namespace Cantera {
m_ScaleSolnNormToResNorm = right.m_ScaleSolnNormToResNorm;
jacCopy_ = right.jacCopy_;
deltax_cp_ = right.deltax_cp_;
deltaX_CP_ = right.deltaX_CP_;
deltaX_Newton_ = right.deltaX_Newton_;
RJd_norm_ = right.RJd_norm_;
lambda_ = right.lambda_;
Jd_ = right.Jd_;
trustDeltaX_ = right.trustDeltaX_;
deltaX_trust_ = right.deltaX_trust_;
trustDelta_ = right.trustDelta_;
Nuu_ = right.Nuu_;
dist_R0_ = right.dist_R0_;
dist_R1_ = right.dist_R1_;
dist_R2_ = right.dist_R2_;
dist_Total_ = right.dist_Total_;
JdJd_norm_ = right.JdJd_norm_;
normTrust_Newton_ = right.normTrust_Newton_;
normTrust_CP_ = right.normTrust_CP_;
return *this;
}
//====================================================================================================================
@ -765,7 +789,7 @@ namespace Cantera {
*
*/
for (int j = 0; j < neq_; j++) {
deltax_cp_[j] = 0.0;
deltaX_CP_[j] = 0.0;
double colFac = 1.0;
if (m_colScaling) {
colFac = 1.0 / m_colScales[j];
@ -774,7 +798,7 @@ namespace Cantera {
if (m_rowScaling) {
rowFac = 1.0 / m_rowScales[i];
}
deltax_cp_[j] -= m_resid[i] * jac.value(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
deltaX_CP_[j] -= m_resid[i] * jac.value(i,j) * colFac * rowFac * m_ewt[j] * m_ewt[j]
/ (m_residWts[i] * m_residWts[i]);
}
}
@ -786,7 +810,7 @@ namespace Cantera {
rowFac = 1.0;
}
for (int j = 0; j < neq_; j++) {
Jd_[i] += deltax_cp_[j] * jac.value(i,j) * rowFac/ m_residWts[i];
Jd_[i] += deltaX_CP_[j] * jac.value(i,j) * rowFac/ m_residWts[i];
}
}
@ -799,7 +823,7 @@ namespace Cantera {
lambda_ = - RJd_norm_ / (JdJd_norm_);
for (int i = 0; i < neq_; i++) {
deltax_cp_[i] *= lambda_;
deltaX_CP_[i] *= lambda_;
}
residNorm2Cauchy_ = m_normResid0 * m_normResid0 - RJd_norm_ * RJd_norm_ / (JdJd_norm_);
@ -814,7 +838,7 @@ namespace Cantera {
}
// Compute the weighted norm of the undamped step size descentDir_[]
normSoln = solnErrorNorm(DATA_PTR(deltax_cp_), "SteepestDescentDir", 10);
normSoln = solnErrorNorm(DATA_PTR(deltaX_CP_), "SteepestDescentDir", 10);
printf("\t\t\tdoCauchyPointSolve: Steepest descent to Cauchy point: \n");
printf("\t\t\t R0 = %g \n", m_normResid0);
@ -832,9 +856,9 @@ namespace Cantera {
int info;
double ff = 1.0E-5;
double *y1 = DATA_PTR(m_wksp);
double s1 = solnErrorNorm(DATA_PTR(deltax_cp_));
double s1 = solnErrorNorm(DATA_PTR(deltaX_CP_));
for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + ff * deltax_cp_[i];
y1[i] = m_y_n[i] + ff * deltaX_CP_[i];
}
/*
* Calculate the residual that would result if y1[] were the new solution vector
@ -940,7 +964,9 @@ namespace Cantera {
/*
* Calculate the trust distances
*/
normTrust_Newton_ = calcTrustDistance(deltaX_Newton_);
normTrust_CP_ = calcTrustDistance(deltaX_CP_);
}
//====================================================================================================================
@ -1040,9 +1066,9 @@ namespace Cantera {
for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) {
double alpha = alphaT[iteration];
for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + alpha * deltax_cp_[i];
y1[i] = m_y_n[i] + alpha * deltaX_CP_[i];
}
sLen = alpha * solnErrorNorm(DATA_PTR(deltax_cp_));
sLen = alpha * solnErrorNorm(DATA_PTR(deltaX_CP_));
/*
* Calculate the residual that would result if y1[] were the new solution vector
* -> m_resid[] contains the result of the residual calculation
@ -1065,7 +1091,7 @@ namespace Cantera {
for (int iteration = 0; iteration < (int) alphaT.size(); iteration++) {
double alpha = alphaT[iteration];
for (int i = 0; i < neq_; i++) {
y1[i] = m_y_n[i] + (1.0 - alpha) * deltax_cp_[i];
y1[i] = m_y_n[i] + (1.0 - alpha) * deltaX_CP_[i];
y1[i] += alpha * Nuu_ * newtDir[i];
}
/*
@ -1270,15 +1296,15 @@ namespace Cantera {
double fabsy;
// we use the old value of the trust region as an indicator
for (int i = 0; i < neq_; i++) {
oldVal = trustDeltaX_[i];
oldVal = deltaX_trust_[i];
fabsy = fabs(m_y_n[i]);
// First off make sure that each trust region vector is 1/2 the size of each variable or smaller
// unless overridden by the deltaStepMininum value.
if (oldVal > 0.5 * fabsy) {
if (fabsy > m_deltaStepMinimum[i]) {
trustDeltaX_[i] = 0.5 * fabsy;
deltaX_trust_[i] = 0.5 * fabsy;
} else {
trustDeltaX_[i] = m_deltaStepMinimum[i];
deltaX_trust_[i] = m_deltaStepMinimum[i];
}
} else {
double newValue = trustDeltaEach * m_ewt[i] / wtSum;
@ -1287,9 +1313,9 @@ namespace Cantera {
} else if (newValue < 0.5 * oldVal) {
newValue = 0.5 * oldVal;
}
trustDeltaX_[i] = newValue;
if (trustDeltaX_[i] > 0.75 * m_deltaStepMaximum[i]) {
trustDeltaX_[i] = m_deltaStepMaximum[i];
deltaX_trust_[i] = newValue;
if (deltaX_trust_[i] > 0.75 * m_deltaStepMaximum[i]) {
deltaX_trust_[i] = m_deltaStepMaximum[i];
}
}
}
@ -1298,10 +1324,10 @@ namespace Cantera {
// Final renormalization.
double sum = 0.0;
for (int i = 0; i < neq_; i++) {
sum += trustDeltaX_[i];
sum += deltaX_trust_[i];
}
for (int i = 0; i < neq_; i++) {
trustDeltaX_[i] = trustDelta_ / sum;
deltaX_trust_[i] = deltaX_trust_[i] / sum;
}
trustDelta_ = 1.0;
@ -1312,12 +1338,51 @@ namespace Cantera {
doublereal sum = 0.0;
doublereal tmp = 0.0;
for (int i = 0; i < neq_; i++) {
tmp = deltaX[i] / trustDeltaX_[i];
tmp = deltaX[i] / deltaX_trust_[i];
sum += tmp * tmp;
}
sum = sqrt(sum / neq_);
return sum;
}
//====================================================================================================================
int NonlinearSolver::calcTrustIntersection(double trustDelta, double &lambda, double &alpha) const
{
double dist;
if (normTrust_Newton_ < trustDelta) {
lambda = 1.0;
alpha = 1.0;
return 2;
}
if (normTrust_Newton_ * Nuu_ > trustDelta) {
alpha = (trustDelta - normTrust_Newton_ * Nuu_) / (normTrust_Newton_ - normTrust_Newton_ * Nuu_);
dist = dist_R0_ + dist_R1_ + alpha * dist_R2_;
lambda = dist / dist_Total_;
return 2;
}
if (normTrust_CP_ > trustDelta) {
lambda = 1.0;
dist = dist_R0_ * trustDelta / normTrust_CP_;
lambda = dist / dist_Total_;
alpha = trustDelta / normTrust_CP_;
return 0;
}
double sumv = 0.0;
for (int i = 0; i < neq_; i++) {
sumv += (deltaX_Newton_[i] / deltaX_trust_[i]) * (deltaX_CP_[i] / deltaX_trust_[i]);
}
double a = normTrust_Newton_ * normTrust_Newton_ * Nuu_ * Nuu_;
double b = 2.0 * Nuu_ * sumv;
double c = normTrust_CP_ * normTrust_CP_ - trustDelta * trustDelta;
alpha =( -b + sqrt( b * b - 4.0 * a * c)) / (2.0 * a);
// alpha = (trustDelta - normTrust_CP_) / (normTrust_Newton_ * Nuu_ - normTrust_CP_);
dist = dist_R0_ + alpha * dist_R1_;
lambda = dist / dist_Total_;
return 1;
}
//====================================================================================================================
/*
*
@ -1625,6 +1690,50 @@ namespace Cantera {
}
return -2;
}
//====================================================================================================================
int NonlinearSolver::dampDogLeg(const doublereal time_curr, const double* y0,
const doublereal *ydot0, const double* step0,
double* const y1, double* const ydot1, double* step1,
double& s1, SquareMatrix& jac, int& loglevel, bool writetitle,
int& num_backtracks)
{
double lambda;
double alpha;
//--------------------------------------------
// Attempt damped step
//--------------------------------------------
// damping coefficient starts at 1.0
m_dampRes = 1.0;
int j, m;
doublereal ff = m_dampBound;
num_backtracks = 0;
/*
* Find the initial value of lambda that satisfies the trust distance
*/
int leg = calcTrustIntersection(trustDelta_, lambda, alpha);
for (m = 0; m < NDAMP; m++) {
}
}
//====================================================================================================================
/*
*

21
Cantera/src/numerics/NonlinearSolver.h
View file

@ -248,6 +248,7 @@ namespace Cantera {
*/
doublereal calcTrustDistance(std::vector<doublereal> const & deltaX) const;
int calcTrustIntersection(double trustDelta, double &lambda, double &alpha) const;
public:
//! Bound the step
/*!
@ -567,6 +568,14 @@ namespace Cantera {
*/
int lambdaToLeg(const double lambda, double &alpha) const;
int calcTrustIntersection(double trustVal, const double &lambda, double &alpha) const;
int dampDogLeg(const doublereal time_curr, const double* y0,
const doublereal *ydot0, const double* step0,
double* const y1, double* const ydot1, double* step1,
double& s1, SquareMatrix& jac, int& loglevel, bool writetitle,
int& num_backtracks);
//! Calculated the expected residual along the double dogleg curve.
/*!
* @param leg 0, 1, or 2 representing the curves of the dogleg
@ -801,7 +810,7 @@ namespace Cantera {
*********************************************************************************************/
//! Steepest descent direction. This is also the distance to the Cauchy Point
std::vector<doublereal> deltax_cp_;
std::vector<doublereal> deltaX_CP_;
//! Newton Step - This is the newton step determined from the straight Jacobian
/*
@ -822,7 +831,7 @@ namespace Cantera {
std::vector<doublereal> Jd_;
//! Vector of trust region values.
std::vector<doublereal> trustDeltaX_;
std::vector<doublereal> deltaX_trust_;
//! Current value of trust radius. This is used with trustDeltaX_ to
//! calculate the max step size.
@ -837,6 +846,14 @@ namespace Cantera {
doublereal dist_Total_;
doublereal JdJd_norm_;
//! Norm of the Newton Step wrt trust region
doublereal normTrust_Newton_;
//! Norm of the Cauchy Step direction wrt trust region
doublereal normTrust_CP_;
/*******************************************************************************************
* OTHER COUNTERS
*****************************************************************************************/

4
ext/f2c_lapack/Makefile.in
View file

@ -100,6 +100,10 @@ dormqr.o \
drscl.o \
dtrtri.o \
dtrti2.o \
dtrtrs.o \
dgecon.o \
dgeequ.o \
dgerfs.o \
ieeeck.o \
ilaenv.o

228
ext/f2c_lapack/dgecon.c Normal file
View file

@ -0,0 +1,228 @@
/* dgecon.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static integer c__1 = 1;
/* Subroutine */ int dgecon_(char *norm, integer *n, doublereal *a, integer *
lda, doublereal *anorm, doublereal *rcond, doublereal *work, integer *
iwork, integer *info, ftnlen norm_len)
{
/* System generated locals */
integer a_dim1, a_offset, i__1;
doublereal d__1;
/* Local variables */
static doublereal sl;
static integer ix;
static doublereal su;
static integer kase, kase1;
static doublereal scale;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int drscl_(integer *, doublereal *, doublereal *,
integer *);
extern doublereal dlamch_(char *, ftnlen);
extern /* Subroutine */ int dlacon_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
extern integer idamax_(integer *, doublereal *, integer *);
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
static doublereal ainvnm;
extern /* Subroutine */ int dlatrs_(char *, char *, char *, char *,
integer *, doublereal *, integer *, doublereal *, doublereal *,
doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen);
static logical onenrm;
static char normin[1];
static doublereal smlnum;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* February 29, 1992 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGECON estimates the reciprocal of the condition number of a general */
/* real matrix A, in either the 1-norm or the infinity-norm, using */
/* the LU factorization computed by DGETRF. */
/* An estimate is obtained for norm(inv(A)), and the reciprocal of the */
/* condition number is computed as */
/* RCOND = 1 / ( norm(A) * norm(inv(A)) ). */
/* Arguments */
/* ========= */
/* NORM (input) CHARACTER*1 */
/* Specifies whether the 1-norm condition number or the */
/* infinity-norm condition number is required: */
/* = '1' or 'O': 1-norm; */
/* = 'I': Infinity-norm. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The factors L and U from the factorization A = P*L*U */
/* as computed by DGETRF. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* ANORM (input) DOUBLE PRECISION */
/* If NORM = '1' or 'O', the 1-norm of the original matrix A. */
/* If NORM = 'I', the infinity-norm of the original matrix A. */
/* RCOND (output) DOUBLE PRECISION */
/* The reciprocal of the condition number of the matrix A, */
/* computed as RCOND = 1/(norm(A) * norm(inv(A))). */
/* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) */
/* IWORK (workspace) INTEGER array, dimension (N) */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--work;
--iwork;
/* Function Body */
*info = 0;
onenrm = *(unsigned char *)norm == '1' || lsame_(norm, "O", (ftnlen)1, (
ftnlen)1);
if (! onenrm && ! lsame_(norm, "I", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*lda < max(1,*n)) {
*info = -4;
} else if (*anorm < 0.) {
*info = -5;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DGECON", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
*rcond = 0.;
if (*n == 0) {
*rcond = 1.;
return 0;
} else if (*anorm == 0.) {
return 0;
}
smlnum = dlamch_("Safe minimum", (ftnlen)12);
/* Estimate the norm of inv(A). */
ainvnm = 0.;
*(unsigned char *)normin = 'N';
if (onenrm) {
kase1 = 1;
} else {
kase1 = 2;
}
kase = 0;
L10:
dlacon_(n, &work[*n + 1], &work[1], &iwork[1], &ainvnm, &kase);
if (kase != 0) {
if (kase == kase1) {
/* Multiply by inv(L). */
dlatrs_("Lower", "No transpose", "Unit", normin, n, &a[a_offset],
lda, &work[1], &sl, &work[(*n << 1) + 1], info, (ftnlen)5,
(ftnlen)12, (ftnlen)4, (ftnlen)1);
/* Multiply by inv(U). */
dlatrs_("Upper", "No transpose", "Non-unit", normin, n, &a[
a_offset], lda, &work[1], &su, &work[*n * 3 + 1], info, (
ftnlen)5, (ftnlen)12, (ftnlen)8, (ftnlen)1);
} else {
/* Multiply by inv(U'). */
dlatrs_("Upper", "Transpose", "Non-unit", normin, n, &a[a_offset],
lda, &work[1], &su, &work[*n * 3 + 1], info, (ftnlen)5, (
ftnlen)9, (ftnlen)8, (ftnlen)1);
/* Multiply by inv(L'). */
dlatrs_("Lower", "Transpose", "Unit", normin, n, &a[a_offset],
lda, &work[1], &sl, &work[(*n << 1) + 1], info, (ftnlen)5,
(ftnlen)9, (ftnlen)4, (ftnlen)1);
}
/* Divide X by 1/(SL*SU) if doing so will not cause overflow. */
scale = sl * su;
*(unsigned char *)normin = 'Y';
if (scale != 1.) {
ix = idamax_(n, &work[1], &c__1);
if (scale < (d__1 = work[ix], abs(d__1)) * smlnum || scale == 0.)
{
goto L20;
}
drscl_(n, &scale, &work[1], &c__1);
}
goto L10;
}
/* Compute the estimate of the reciprocal condition number. */
if (ainvnm != 0.) {
*rcond = 1. / ainvnm / *anorm;
}
L20:
return 0;
/* End of DGECON */
} /* dgecon_ */

297
ext/f2c_lapack/dgeequ.c Normal file
View file

@ -0,0 +1,297 @@
/* dgeequ.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Subroutine */ int dgeequ_(integer *m, integer *n, doublereal *a, integer *
lda, doublereal *r__, doublereal *c__, doublereal *rowcnd, doublereal
*colcnd, doublereal *amax, integer *info)
{
/* System generated locals */
integer a_dim1, a_offset, i__1, i__2;
doublereal d__1, d__2, d__3;
/* Local variables */
static integer i__, j;
static doublereal rcmin, rcmax;
extern doublereal dlamch_(char *, ftnlen);
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen);
static doublereal bignum, smlnum;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* March 31, 1993 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGEEQU computes row and column scalings intended to equilibrate an */
/* M-by-N matrix A and reduce its condition number. R returns the row */
/* scale factors and C the column scale factors, chosen to try to make */
/* the largest element in each row and column of the matrix B with */
/* elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1. */
/* R(i) and C(j) are restricted to be between SMLNUM = smallest safe */
/* number and BIGNUM = largest safe number. Use of these scaling */
/* factors is not guaranteed to reduce the condition number of A but */
/* works well in practice. */
/* Arguments */
/* ========= */
/* M (input) INTEGER */
/* The number of rows of the matrix A. M >= 0. */
/* N (input) INTEGER */
/* The number of columns of the matrix A. N >= 0. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The M-by-N matrix whose equilibration factors are */
/* to be computed. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,M). */
/* R (output) DOUBLE PRECISION array, dimension (M) */
/* If INFO = 0 or INFO > M, R contains the row scale factors */
/* for A. */
/* C (output) DOUBLE PRECISION array, dimension (N) */
/* If INFO = 0, C contains the column scale factors for A. */
/* ROWCND (output) DOUBLE PRECISION */
/* If INFO = 0 or INFO > M, ROWCND contains the ratio of the */
/* smallest R(i) to the largest R(i). If ROWCND >= 0.1 and */
/* AMAX is neither too large nor too small, it is not worth */
/* scaling by R. */
/* COLCND (output) DOUBLE PRECISION */
/* If INFO = 0, COLCND contains the ratio of the smallest */
/* C(i) to the largest C(i). If COLCND >= 0.1, it is not */
/* worth scaling by C. */
/* AMAX (output) DOUBLE PRECISION */
/* Absolute value of largest matrix element. If AMAX is very */
/* close to overflow or very close to underflow, the matrix */
/* should be scaled. */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: if INFO = i, and i is */
/* <= M: the i-th row of A is exactly zero */
/* > M: the (i-M)-th column of A is exactly zero */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
--r__;
--c__;
/* Function Body */
*info = 0;
if (*m < 0) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*lda < max(1,*m)) {
*info = -4;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DGEEQU", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
if (*m == 0 || *n == 0) {
*rowcnd = 1.;
*colcnd = 1.;
*amax = 0.;
return 0;
}
/* Get machine constants. */
smlnum = dlamch_("S", (ftnlen)1);
bignum = 1. / smlnum;
/* Compute row scale factors. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
r__[i__] = 0.;
/* L10: */
}
/* Find the maximum element in each row. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
d__2 = r__[i__], d__3 = (d__1 = a[i__ + j * a_dim1], abs(d__1));
r__[i__] = max(d__2,d__3);
/* L20: */
}
/* L30: */
}
/* Find the maximum and minimum scale factors. */
rcmin = bignum;
rcmax = 0.;
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MAX */
d__1 = rcmax, d__2 = r__[i__];
rcmax = max(d__1,d__2);
/* Computing MIN */
d__1 = rcmin, d__2 = r__[i__];
rcmin = min(d__1,d__2);
/* L40: */
}
*amax = rcmax;
if (rcmin == 0.) {
/* Find the first zero scale factor and return an error code. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
if (r__[i__] == 0.) {
*info = i__;
return 0;
}
/* L50: */
}
} else {
/* Invert the scale factors. */
i__1 = *m;
for (i__ = 1; i__ <= i__1; ++i__) {
/* Computing MIN */
/* Computing MAX */
d__2 = r__[i__];
d__1 = max(d__2,smlnum);
r__[i__] = 1. / min(d__1,bignum);
/* L60: */
}
/* Compute ROWCND = min(R(I)) / max(R(I)) */
*rowcnd = max(rcmin,smlnum) / min(rcmax,bignum);
}
/* Compute column scale factors */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
c__[j] = 0.;
/* L70: */
}
/* Find the maximum element in each column, */
/* assuming the row scaling computed above. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
i__2 = *m;
for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
d__2 = c__[j], d__3 = (d__1 = a[i__ + j * a_dim1], abs(d__1)) *
r__[i__];
c__[j] = max(d__2,d__3);
/* L80: */
}
/* L90: */
}
/* Find the maximum and minimum scale factors. */
rcmin = bignum;
rcmax = 0.;
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
d__1 = rcmin, d__2 = c__[j];
rcmin = min(d__1,d__2);
/* Computing MAX */
d__1 = rcmax, d__2 = c__[j];
rcmax = max(d__1,d__2);
/* L100: */
}
if (rcmin == 0.) {
/* Find the first zero scale factor and return an error code. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
if (c__[j] == 0.) {
*info = *m + j;
return 0;
}
/* L110: */
}
} else {
/* Invert the scale factors. */
i__1 = *n;
for (j = 1; j <= i__1; ++j) {
/* Computing MIN */
/* Computing MAX */
d__2 = c__[j];
d__1 = max(d__2,smlnum);
c__[j] = 1. / min(d__1,bignum);
/* L120: */
}
/* Compute COLCND = min(C(J)) / max(C(J)) */
*colcnd = max(rcmin,smlnum) / min(rcmax,bignum);
}
return 0;
/* End of DGEEQU */
} /* dgeequ_ */

420
ext/f2c_lapack/dgerfs.c Normal file
View file

@ -0,0 +1,420 @@
/* dgerfs.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static integer c__1 = 1;
static doublereal c_b15 = -1.;
static doublereal c_b17 = 1.;
/* Subroutine */ int dgerfs_(char *trans, integer *n, integer *nrhs,
doublereal *a, integer *lda, doublereal *af, integer *ldaf, integer *
ipiv, doublereal *b, integer *ldb, doublereal *x, integer *ldx,
doublereal *ferr, doublereal *berr, doublereal *work, integer *iwork,
integer *info, ftnlen trans_len)
{
/* System generated locals */
integer a_dim1, a_offset, af_dim1, af_offset, b_dim1, b_offset, x_dim1,
x_offset, i__1, i__2, i__3;
doublereal d__1, d__2, d__3;
/* Local variables */
static integer i__, j, k;
static doublereal s, xk;
static integer nz;
static doublereal eps;
static integer kase;
static doublereal safe1, safe2;
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int dgemv_(char *, integer *, integer *,
doublereal *, doublereal *, integer *, doublereal *, integer *,
doublereal *, doublereal *, integer *, ftnlen), dcopy_(integer *,
doublereal *, integer *, doublereal *, integer *), daxpy_(integer
*, doublereal *, doublereal *, integer *, doublereal *, integer *)
;
static integer count;
extern doublereal dlamch_(char *, ftnlen);
extern /* Subroutine */ int dlacon_(integer *, doublereal *, doublereal *,
integer *, doublereal *, integer *);
static doublereal safmin;
extern /* Subroutine */ int xerbla_(char *, integer *, ftnlen), dgetrs_(
char *, integer *, integer *, doublereal *, integer *, integer *,
doublereal *, integer *, integer *, ftnlen);
static logical notran;
static char transt[1];
static doublereal lstres;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* September 30, 1994 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DGERFS improves the computed solution to a system of linear */
/* equations and provides error bounds and backward error estimates for */
/* the solution. */
/* Arguments */
/* ========= */
/* TRANS (input) CHARACTER*1 */
/* Specifies the form of the system of equations: */
/* = 'N': A * X = B (No transpose) */
/* = 'T': A**T * X = B (Transpose) */
/* = 'C': A**H * X = B (Conjugate transpose = Transpose) */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* NRHS (input) INTEGER */
/* The number of right hand sides, i.e., the number of columns */
/* of the matrices B and X. NRHS >= 0. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The original N-by-N matrix A. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* AF (input) DOUBLE PRECISION array, dimension (LDAF,N) */
/* The factors L and U from the factorization A = P*L*U */
/* as computed by DGETRF. */
/* LDAF (input) INTEGER */
/* The leading dimension of the array AF. LDAF >= max(1,N). */
/* IPIV (input) INTEGER array, dimension (N) */
/* The pivot indices from DGETRF; for 1<=i<=N, row i of the */
/* matrix was interchanged with row IPIV(i). */
/* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS) */
/* The right hand side matrix B. */
/* LDB (input) INTEGER */
/* The leading dimension of the array B. LDB >= max(1,N). */
/* X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS) */
/* On entry, the solution matrix X, as computed by DGETRS. */
/* On exit, the improved solution matrix X. */
/* LDX (input) INTEGER */
/* The leading dimension of the array X. LDX >= max(1,N). */
/* FERR (output) DOUBLE PRECISION array, dimension (NRHS) */
/* The estimated forward error bound for each solution vector */
/* X(j) (the j-th column of the solution matrix X). */
/* If XTRUE is the true solution corresponding to X(j), FERR(j) */
/* is an estimated upper bound for the magnitude of the largest */
/* element in (X(j) - XTRUE) divided by the magnitude of the */
/* largest element in X(j). The estimate is as reliable as */
/* the estimate for RCOND, and is almost always a slight */
/* overestimate of the true error. */
/* BERR (output) DOUBLE PRECISION array, dimension (NRHS) */
/* The componentwise relative backward error of each solution */
/* vector X(j) (i.e., the smallest relative change in */
/* any element of A or B that makes X(j) an exact solution). */
/* WORK (workspace) DOUBLE PRECISION array, dimension (3*N) */
/* IWORK (workspace) INTEGER array, dimension (N) */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* Internal Parameters */
/* =================== */
/* ITMAX is the maximum number of steps of iterative refinement. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
af_dim1 = *ldaf;
af_offset = 1 + af_dim1;
af -= af_offset;
--ipiv;
b_dim1 = *ldb;
b_offset = 1 + b_dim1;
b -= b_offset;
x_dim1 = *ldx;
x_offset = 1 + x_dim1;
x -= x_offset;
--ferr;
--berr;
--work;
--iwork;
/* Function Body */
*info = 0;
notran = lsame_(trans, "N", (ftnlen)1, (ftnlen)1);
if (! notran && ! lsame_(trans, "T", (ftnlen)1, (ftnlen)1) && ! lsame_(
trans, "C", (ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (*n < 0) {
*info = -2;
} else if (*nrhs < 0) {
*info = -3;
} else if (*lda < max(1,*n)) {
*info = -5;
} else if (*ldaf < max(1,*n)) {
*info = -7;
} else if (*ldb < max(1,*n)) {
*info = -10;
} else if (*ldx < max(1,*n)) {
*info = -12;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DGERFS", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
if (*n == 0 || *nrhs == 0) {
i__1 = *nrhs;
for (j = 1; j <= i__1; ++j) {
ferr[j] = 0.;
berr[j] = 0.;
/* L10: */
}
return 0;
}
if (notran) {
*(unsigned char *)transt = 'T';
} else {
*(unsigned char *)transt = 'N';
}
/* NZ = maximum number of nonzero elements in each row of A, plus 1 */
nz = *n + 1;
eps = dlamch_("Epsilon", (ftnlen)7);
safmin = dlamch_("Safe minimum", (ftnlen)12);
safe1 = nz * safmin;
safe2 = safe1 / eps;
/* Do for each right hand side */
i__1 = *nrhs;
for (j = 1; j <= i__1; ++j) {
count = 1;
lstres = 3.;
L20:
/* Loop until stopping criterion is satisfied. */
/* Compute residual R = B - op(A) * X, */
/* where op(A) = A, A**T, or A**H, depending on TRANS. */
dcopy_(n, &b[j * b_dim1 + 1], &c__1, &work[*n + 1], &c__1);
dgemv_(trans, n, n, &c_b15, &a[a_offset], lda, &x[j * x_dim1 + 1], &
c__1, &c_b17, &work[*n + 1], &c__1, (ftnlen)1);
/* Compute componentwise relative backward error from formula */
/* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) ) */
/* where abs(Z) is the componentwise absolute value of the matrix */
/* or vector Z. If the i-th component of the denominator is less */
/* than SAFE2, then SAFE1 is added to the i-th components of the */
/* numerator and denominator before dividing. */
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
work[i__] = (d__1 = b[i__ + j * b_dim1], abs(d__1));
/* L30: */
}
/* Compute abs(op(A))*abs(X) + abs(B). */
if (notran) {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
xk = (d__1 = x[k + j * x_dim1], abs(d__1));
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
work[i__] += (d__1 = a[i__ + k * a_dim1], abs(d__1)) * xk;
/* L40: */
}
/* L50: */
}
} else {
i__2 = *n;
for (k = 1; k <= i__2; ++k) {
s = 0.;
i__3 = *n;
for (i__ = 1; i__ <= i__3; ++i__) {
s += (d__1 = a[i__ + k * a_dim1], abs(d__1)) * (d__2 = x[
i__ + j * x_dim1], abs(d__2));
/* L60: */
}
work[k] += s;
/* L70: */
}
}
s = 0.;
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
if (work[i__] > safe2) {
/* Computing MAX */
d__2 = s, d__3 = (d__1 = work[*n + i__], abs(d__1)) / work[
i__];
s = max(d__2,d__3);
} else {
/* Computing MAX */
d__2 = s, d__3 = ((d__1 = work[*n + i__], abs(d__1)) + safe1)
/ (work[i__] + safe1);
s = max(d__2,d__3);
}
/* L80: */
}
berr[j] = s;
/* Test stopping criterion. Continue iterating if */
/* 1) The residual BERR(J) is larger than machine epsilon, and */
/* 2) BERR(J) decreased by at least a factor of 2 during the */
/* last iteration, and */
/* 3) At most ITMAX iterations tried. */
if (berr[j] > eps && berr[j] * 2. <= lstres && count <= 5) {
/* Update solution and try again. */
dgetrs_(trans, n, &c__1, &af[af_offset], ldaf, &ipiv[1], &work[*n
+ 1], n, info, (ftnlen)1);
daxpy_(n, &c_b17, &work[*n + 1], &c__1, &x[j * x_dim1 + 1], &c__1)
;
lstres = berr[j];
++count;
goto L20;
}
/* Bound error from formula */
/* norm(X - XTRUE) / norm(X) .le. FERR = */
/* norm( abs(inv(op(A)))* */
/* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X) */
/* where */
/* norm(Z) is the magnitude of the largest component of Z */
/* inv(op(A)) is the inverse of op(A) */
/* abs(Z) is the componentwise absolute value of the matrix or */
/* vector Z */
/* NZ is the maximum number of nonzeros in any row of A, plus 1 */
/* EPS is machine epsilon */
/* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B)) */
/* is incremented by SAFE1 if the i-th component of */
/* abs(op(A))*abs(X) + abs(B) is less than SAFE2. */
/* Use DLACON to estimate the infinity-norm of the matrix */
/* inv(op(A)) * diag(W), */
/* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) */
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
if (work[i__] > safe2) {
work[i__] = (d__1 = work[*n + i__], abs(d__1)) + nz * eps *
work[i__];
} else {
work[i__] = (d__1 = work[*n + i__], abs(d__1)) + nz * eps *
work[i__] + safe1;
}
/* L90: */
}
kase = 0;
L100:
dlacon_(n, &work[(*n << 1) + 1], &work[*n + 1], &iwork[1], &ferr[j], &
kase);
if (kase != 0) {
if (kase == 1) {
/* Multiply by diag(W)*inv(op(A)**T). */
dgetrs_(transt, n, &c__1, &af[af_offset], ldaf, &ipiv[1], &
work[*n + 1], n, info, (ftnlen)1);
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
work[*n + i__] = work[i__] * work[*n + i__];
/* L110: */
}
} else {
/* Multiply by inv(op(A))*diag(W). */
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
work[*n + i__] = work[i__] * work[*n + i__];
/* L120: */
}
dgetrs_(trans, n, &c__1, &af[af_offset], ldaf, &ipiv[1], &
work[*n + 1], n, info, (ftnlen)1);
}
goto L100;
}
/* Normalize error. */
lstres = 0.;
i__2 = *n;
for (i__ = 1; i__ <= i__2; ++i__) {
/* Computing MAX */
d__2 = lstres, d__3 = (d__1 = x[i__ + j * x_dim1], abs(d__1));
lstres = max(d__2,d__3);
/* L130: */
}
if (lstres != 0.) {
ferr[j] /= lstres;
}
/* L140: */
}
return 0;
/* End of DGERFS */
} /* dgerfs_ */

187
ext/f2c_lapack/dtrtrs.c Normal file
View file

@ -0,0 +1,187 @@
/* dtrtrs.f -- translated by f2c (version 20031025).
You must link the resulting object file with libf2c:
on Microsoft Windows system, link with libf2c.lib;
on Linux or Unix systems, link with .../path/to/libf2c.a -lm
or, if you install libf2c.a in a standard place, with -lf2c -lm
-- in that order, at the end of the command line, as in
cc *.o -lf2c -lm
Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
http://www.netlib.org/f2c/libf2c.zip
*/
#include "f2c.h"
/* Table of constant values */
static doublereal c_b12 = 1.;
/* Subroutine */ int dtrtrs_(char *uplo, char *trans, char *diag, integer *n,
integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *
ldb, integer *info, ftnlen uplo_len, ftnlen trans_len, ftnlen
diag_len)
{
/* System generated locals */
integer a_dim1, a_offset, b_dim1, b_offset, i__1;
/* Local variables */
extern logical lsame_(char *, char *, ftnlen, ftnlen);
extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
integer *, integer *, doublereal *, doublereal *, integer *,
doublereal *, integer *, ftnlen, ftnlen, ftnlen, ftnlen), xerbla_(
char *, integer *, ftnlen);
static logical nounit;
/* -- LAPACK routine (version 3.0) -- */
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
/* Courant Institute, Argonne National Lab, and Rice University */
/* March 31, 1993 */
/* .. Scalar Arguments .. */
/* .. */
/* .. Array Arguments .. */
/* .. */
/* Purpose */
/* ======= */
/* DTRTRS solves a triangular system of the form */
/* A * X = B or A**T * X = B, */
/* where A is a triangular matrix of order N, and B is an N-by-NRHS */
/* matrix. A check is made to verify that A is nonsingular. */
/* Arguments */
/* ========= */
/* UPLO (input) CHARACTER*1 */
/* = 'U': A is upper triangular; */
/* = 'L': A is lower triangular. */
/* TRANS (input) CHARACTER*1 */
/* Specifies the form of the system of equations: */
/* = 'N': A * X = B (No transpose) */
/* = 'T': A**T * X = B (Transpose) */
/* = 'C': A**H * X = B (Conjugate transpose = Transpose) */
/* DIAG (input) CHARACTER*1 */
/* = 'N': A is non-unit triangular; */
/* = 'U': A is unit triangular. */
/* N (input) INTEGER */
/* The order of the matrix A. N >= 0. */
/* NRHS (input) INTEGER */
/* The number of right hand sides, i.e., the number of columns */
/* of the matrix B. NRHS >= 0. */
/* A (input) DOUBLE PRECISION array, dimension (LDA,N) */
/* The triangular matrix A. If UPLO = 'U', the leading N-by-N */
/* upper triangular part of the array A contains the upper */
/* triangular matrix, and the strictly lower triangular part of */
/* A is not referenced. If UPLO = 'L', the leading N-by-N lower */
/* triangular part of the array A contains the lower triangular */
/* matrix, and the strictly upper triangular part of A is not */
/* referenced. If DIAG = 'U', the diagonal elements of A are */
/* also not referenced and are assumed to be 1. */
/* LDA (input) INTEGER */
/* The leading dimension of the array A. LDA >= max(1,N). */
/* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) */
/* On entry, the right hand side matrix B. */
/* On exit, if INFO = 0, the solution matrix X. */
/* LDB (input) INTEGER */
/* The leading dimension of the array B. LDB >= max(1,N). */
/* INFO (output) INTEGER */
/* = 0: successful exit */
/* < 0: if INFO = -i, the i-th argument had an illegal value */
/* > 0: if INFO = i, the i-th diagonal element of A is zero, */
/* indicating that the matrix is singular and the solutions */
/* X have not been computed. */
/* ===================================================================== */
/* .. Parameters .. */
/* .. */
/* .. Local Scalars .. */
/* .. */
/* .. External Functions .. */
/* .. */
/* .. External Subroutines .. */
/* .. */
/* .. Intrinsic Functions .. */
/* .. */
/* .. Executable Statements .. */
/* Test the input parameters. */
/* Parameter adjustments */
a_dim1 = *lda;
a_offset = 1 + a_dim1;
a -= a_offset;
b_dim1 = *ldb;
b_offset = 1 + b_dim1;
b -= b_offset;
/* Function Body */
*info = 0;
nounit = lsame_(diag, "N", (ftnlen)1, (ftnlen)1);
if (! lsame_(uplo, "U", (ftnlen)1, (ftnlen)1) && ! lsame_(uplo, "L", (
ftnlen)1, (ftnlen)1)) {
*info = -1;
} else if (! lsame_(trans, "N", (ftnlen)1, (ftnlen)1) && ! lsame_(trans,
"T", (ftnlen)1, (ftnlen)1) && ! lsame_(trans, "C", (ftnlen)1, (
ftnlen)1)) {
*info = -2;
} else if (! nounit && ! lsame_(diag, "U", (ftnlen)1, (ftnlen)1)) {
*info = -3;
} else if (*n < 0) {
*info = -4;
} else if (*nrhs < 0) {
*info = -5;
} else if (*lda < max(1,*n)) {
*info = -7;
} else if (*ldb < max(1,*n)) {
*info = -9;
}
if (*info != 0) {
i__1 = -(*info);
xerbla_("DTRTRS", &i__1, (ftnlen)6);
return 0;
}
/* Quick return if possible */
if (*n == 0) {
return 0;
}
/* Check for singularity. */
if (nounit) {
i__1 = *n;
for (*info = 1; *info <= i__1; ++(*info)) {
if (a[*info + *info * a_dim1] == 0.) {
return 0;
}
/* L10: */
}
}
*info = 0;
/* Solve A * x = b or A' * x = b. */
dtrsm_("Left", uplo, trans, diag, n, nrhs, &c_b12, &a[a_offset], lda, &b[
b_offset], ldb, (ftnlen)4, (ftnlen)1, (ftnlen)1, (ftnlen)1);
return 0;
/* End of DTRTRS */
} /* dtrtrs_ */

4
ext/lapack/Makefile.in
View file

@ -67,6 +67,10 @@ dorml2.o \
dormlq.o \
dormqr.o \
drscl.o \
dtrtrs.o \
dgerfs.o \
dgecon.o \
dgeequ.o \
ilaenv.o
#SRCS = $(OBJS:.o=.cpp)

181
ext/lapack/dgecon.f Normal file
View file

@ -0,0 +1,181 @@
SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
$ INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. Scalar Arguments ..
CHARACTER NORM
INTEGER INFO, LDA, N
DOUBLE PRECISION ANORM, RCOND
* ..
* .. Array Arguments ..
INTEGER IWORK( * )
DOUBLE PRECISION A( LDA, * ), WORK( * )
* ..
*
* Purpose
* =======
*
* DGECON estimates the reciprocal of the condition number of a general
* real matrix A, in either the 1-norm or the infinity-norm, using
* the LU factorization computed by DGETRF.
*
* An estimate is obtained for norm(inv(A)), and the reciprocal of the
* condition number is computed as
* RCOND = 1 / ( norm(A) * norm(inv(A)) ).
*
* Arguments
* =========
*
* NORM (input) CHARACTER*1
* Specifies whether the 1-norm condition number or the
* infinity-norm condition number is required:
* = '1' or 'O': 1-norm;
* = 'I': Infinity-norm.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
* The factors L and U from the factorization A = P*L*U
* as computed by DGETRF.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* ANORM (input) DOUBLE PRECISION
* If NORM = '1' or 'O', the 1-norm of the original matrix A.
* If NORM = 'I', the infinity-norm of the original matrix A.
*
* RCOND (output) DOUBLE PRECISION
* The reciprocal of the condition number of the matrix A,
* computed as RCOND = 1/(norm(A) * norm(inv(A))).
*
* WORK (workspace) DOUBLE PRECISION array, dimension (4*N)
*
* IWORK (workspace) INTEGER array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL ONENRM
CHARACTER NORMIN
INTEGER IX, KASE, KASE1
DOUBLE PRECISION AINVNM, SCALE, SL, SMLNUM, SU
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DLAMCH
EXTERNAL LSAME, IDAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DLACON, DLATRS, DRSCL, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -4
ELSE IF( ANORM.LT.ZERO ) THEN
INFO = -5
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGECON', -INFO )
RETURN
END IF
*
* Quick return if possible
*
RCOND = ZERO
IF( N.EQ.0 ) THEN
RCOND = ONE
RETURN
ELSE IF( ANORM.EQ.ZERO ) THEN
RETURN
END IF
*
SMLNUM = DLAMCH( 'Safe minimum' )
*
* Estimate the norm of inv(A).
*
AINVNM = ZERO
NORMIN = 'N'
IF( ONENRM ) THEN
KASE1 = 1
ELSE
KASE1 = 2
END IF
KASE = 0
10 CONTINUE
CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE )
IF( KASE.NE.0 ) THEN
IF( KASE.EQ.KASE1 ) THEN
*
* Multiply by inv(L).
*
CALL DLATRS( 'Lower', 'No transpose', 'Unit', NORMIN, N, A,
$ LDA, WORK, SL, WORK( 2*N+1 ), INFO )
*
* Multiply by inv(U).
*
CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
$ A, LDA, WORK, SU, WORK( 3*N+1 ), INFO )
ELSE
*
* Multiply by inv(U').
*
CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
$ LDA, WORK, SU, WORK( 3*N+1 ), INFO )
*
* Multiply by inv(L').
*
CALL DLATRS( 'Lower', 'Transpose', 'Unit', NORMIN, N, A,
$ LDA, WORK, SL, WORK( 2*N+1 ), INFO )
END IF
*
* Divide X by 1/(SL*SU) if doing so will not cause overflow.
*
SCALE = SL*SU
NORMIN = 'Y'
IF( SCALE.NE.ONE ) THEN
IX = IDAMAX( N, WORK, 1 )
IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
$ GO TO 20
CALL DRSCL( N, SCALE, WORK, 1 )
END IF
GO TO 10
END IF
*
* Compute the estimate of the reciprocal condition number.
*
IF( AINVNM.NE.ZERO )
$ RCOND = ( ONE / AINVNM ) / ANORM
*
20 CONTINUE
RETURN
*
* End of DGECON
*
END

226
ext/lapack/dgeequ.f Normal file
View file

@ -0,0 +1,226 @@
SUBROUTINE DGEEQU( M, N, A, LDA, R, C, ROWCND, COLCND, AMAX,
$ INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* March 31, 1993
*
* .. Scalar Arguments ..
INTEGER INFO, LDA, M, N
DOUBLE PRECISION AMAX, COLCND, ROWCND
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), C( * ), R( * )
* ..
*
* Purpose
* =======
*
* DGEEQU computes row and column scalings intended to equilibrate an
* M-by-N matrix A and reduce its condition number. R returns the row
* scale factors and C the column scale factors, chosen to try to make
* the largest element in each row and column of the matrix B with
* elements B(i,j)=R(i)*A(i,j)*C(j) have absolute value 1.
*
* R(i) and C(j) are restricted to be between SMLNUM = smallest safe
* number and BIGNUM = largest safe number. Use of these scaling
* factors is not guaranteed to reduce the condition number of A but
* works well in practice.
*
* Arguments
* =========
*
* M (input) INTEGER
* The number of rows of the matrix A. M >= 0.
*
* N (input) INTEGER
* The number of columns of the matrix A. N >= 0.
*
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
* The M-by-N matrix whose equilibration factors are
* to be computed.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,M).
*
* R (output) DOUBLE PRECISION array, dimension (M)
* If INFO = 0 or INFO > M, R contains the row scale factors
* for A.
*
* C (output) DOUBLE PRECISION array, dimension (N)
* If INFO = 0, C contains the column scale factors for A.
*
* ROWCND (output) DOUBLE PRECISION
* If INFO = 0 or INFO > M, ROWCND contains the ratio of the
* smallest R(i) to the largest R(i). If ROWCND >= 0.1 and
* AMAX is neither too large nor too small, it is not worth
* scaling by R.
*
* COLCND (output) DOUBLE PRECISION
* If INFO = 0, COLCND contains the ratio of the smallest
* C(i) to the largest C(i). If COLCND >= 0.1, it is not
* worth scaling by C.
*
* AMAX (output) DOUBLE PRECISION
* Absolute value of largest matrix element. If AMAX is very
* close to overflow or very close to underflow, the matrix
* should be scaled.
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, and i is
* <= M: the i-th row of A is exactly zero
* > M: the (i-M)-th column of A is exactly zero
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE, ZERO
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
* ..
* .. Local Scalars ..
INTEGER I, J
DOUBLE PRECISION BIGNUM, RCMAX, RCMIN, SMLNUM
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH
EXTERNAL DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
IF( M.LT.0 ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEEQU', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( M.EQ.0 .OR. N.EQ.0 ) THEN
ROWCND = ONE
COLCND = ONE
AMAX = ZERO
RETURN
END IF
*
* Get machine constants.
*
SMLNUM = DLAMCH( 'S' )
BIGNUM = ONE / SMLNUM
*
* Compute row scale factors.
*
DO 10 I = 1, M
R( I ) = ZERO
10 CONTINUE
*
* Find the maximum element in each row.
*
DO 30 J = 1, N
DO 20 I = 1, M
R( I ) = MAX( R( I ), ABS( A( I, J ) ) )
20 CONTINUE
30 CONTINUE
*
* Find the maximum and minimum scale factors.
*
RCMIN = BIGNUM
RCMAX = ZERO
DO 40 I = 1, M
RCMAX = MAX( RCMAX, R( I ) )
RCMIN = MIN( RCMIN, R( I ) )
40 CONTINUE
AMAX = RCMAX
*
IF( RCMIN.EQ.ZERO ) THEN
*
* Find the first zero scale factor and return an error code.
*
DO 50 I = 1, M
IF( R( I ).EQ.ZERO ) THEN
INFO = I
RETURN
END IF
50 CONTINUE
ELSE
*
* Invert the scale factors.
*
DO 60 I = 1, M
R( I ) = ONE / MIN( MAX( R( I ), SMLNUM ), BIGNUM )
60 CONTINUE
*
* Compute ROWCND = min(R(I)) / max(R(I))
*
ROWCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
END IF
*
* Compute column scale factors
*
DO 70 J = 1, N
C( J ) = ZERO
70 CONTINUE
*
* Find the maximum element in each column,
* assuming the row scaling computed above.
*
DO 90 J = 1, N
DO 80 I = 1, M
C( J ) = MAX( C( J ), ABS( A( I, J ) )*R( I ) )
80 CONTINUE
90 CONTINUE
*
* Find the maximum and minimum scale factors.
*
RCMIN = BIGNUM
RCMAX = ZERO
DO 100 J = 1, N
RCMIN = MIN( RCMIN, C( J ) )
RCMAX = MAX( RCMAX, C( J ) )
100 CONTINUE
*
IF( RCMIN.EQ.ZERO ) THEN
*
* Find the first zero scale factor and return an error code.
*
DO 110 J = 1, N
IF( C( J ).EQ.ZERO ) THEN
INFO = M + J
RETURN
END IF
110 CONTINUE
ELSE
*
* Invert the scale factors.
*
DO 120 J = 1, N
C( J ) = ONE / MIN( MAX( C( J ), SMLNUM ), BIGNUM )
120 CONTINUE
*
* Compute COLCND = min(C(J)) / max(C(J))
*
COLCND = MAX( RCMIN, SMLNUM ) / MIN( RCMAX, BIGNUM )
END IF
*
RETURN
*
* End of DGEEQU
*
END

332
ext/lapack/dgerfs.f Normal file
View file

@ -0,0 +1,332 @@
SUBROUTINE DGERFS( TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, B, LDB,
$ X, LDX, FERR, BERR, WORK, IWORK, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
CHARACTER TRANS
INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS
* ..
* .. Array Arguments ..
INTEGER IPIV( * ), IWORK( * )
DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ),
$ BERR( * ), FERR( * ), WORK( * ), X( LDX, * )
* ..
*
* Purpose
* =======
*
* DGERFS improves the computed solution to a system of linear
* equations and provides error bounds and backward error estimates for
* the solution.
*
* Arguments
* =========
*
* TRANS (input) CHARACTER*1
* Specifies the form of the system of equations:
* = 'N': A * X = B (No transpose)
* = 'T': A**T * X = B (Transpose)
* = 'C': A**H * X = B (Conjugate transpose = Transpose)
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrices B and X. NRHS >= 0.
*
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
* The original N-by-N matrix A.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* AF (input) DOUBLE PRECISION array, dimension (LDAF,N)
* The factors L and U from the factorization A = P*L*U
* as computed by DGETRF.
*
* LDAF (input) INTEGER
* The leading dimension of the array AF. LDAF >= max(1,N).
*
* IPIV (input) INTEGER array, dimension (N)
* The pivot indices from DGETRF; for 1<=i<=N, row i of the
* matrix was interchanged with row IPIV(i).
*
* B (input) DOUBLE PRECISION array, dimension (LDB,NRHS)
* The right hand side matrix B.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* X (input/output) DOUBLE PRECISION array, dimension (LDX,NRHS)
* On entry, the solution matrix X, as computed by DGETRS.
* On exit, the improved solution matrix X.
*
* LDX (input) INTEGER
* The leading dimension of the array X. LDX >= max(1,N).
*
* FERR (output) DOUBLE PRECISION array, dimension (NRHS)
* The estimated forward error bound for each solution vector
* X(j) (the j-th column of the solution matrix X).
* If XTRUE is the true solution corresponding to X(j), FERR(j)
* is an estimated upper bound for the magnitude of the largest
* element in (X(j) - XTRUE) divided by the magnitude of the
* largest element in X(j). The estimate is as reliable as
* the estimate for RCOND, and is almost always a slight
* overestimate of the true error.
*
* BERR (output) DOUBLE PRECISION array, dimension (NRHS)
* The componentwise relative backward error of each solution
* vector X(j) (i.e., the smallest relative change in
* any element of A or B that makes X(j) an exact solution).
*
* WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
*
* IWORK (workspace) INTEGER array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
*
* Internal Parameters
* ===================
*
* ITMAX is the maximum number of steps of iterative refinement.
*
* =====================================================================
*
* .. Parameters ..
INTEGER ITMAX
PARAMETER ( ITMAX = 5 )
DOUBLE PRECISION ZERO
PARAMETER ( ZERO = 0.0D+0 )
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
DOUBLE PRECISION TWO
PARAMETER ( TWO = 2.0D+0 )
DOUBLE PRECISION THREE
PARAMETER ( THREE = 3.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL NOTRAN
CHARACTER TRANST
INTEGER COUNT, I, J, K, KASE, NZ
DOUBLE PRECISION EPS, LSTRES, S, SAFE1, SAFE2, SAFMIN, XK
* ..
* .. External Subroutines ..
EXTERNAL DAXPY, DCOPY, DGEMV, DGETRS, DLACON, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX
* ..
* .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH
EXTERNAL LSAME, DLAMCH
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
NOTRAN = LSAME( TRANS, 'N' )
IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
$ LSAME( TRANS, 'C' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( NRHS.LT.0 ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDAF.LT.MAX( 1, N ) ) THEN
INFO = -7
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -10
ELSE IF( LDX.LT.MAX( 1, N ) ) THEN
INFO = -12
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGERFS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 .OR. NRHS.EQ.0 ) THEN
DO 10 J = 1, NRHS
FERR( J ) = ZERO
BERR( J ) = ZERO
10 CONTINUE
RETURN
END IF
*
IF( NOTRAN ) THEN
TRANST = 'T'
ELSE
TRANST = 'N'
END IF
*
* NZ = maximum number of nonzero elements in each row of A, plus 1
*
NZ = N + 1
EPS = DLAMCH( 'Epsilon' )
SAFMIN = DLAMCH( 'Safe minimum' )
SAFE1 = NZ*SAFMIN
SAFE2 = SAFE1 / EPS
*
* Do for each right hand side
*
DO 140 J = 1, NRHS
*
COUNT = 1
LSTRES = THREE
20 CONTINUE
*
* Loop until stopping criterion is satisfied.
*
* Compute residual R = B - op(A) * X,
* where op(A) = A, A**T, or A**H, depending on TRANS.
*
CALL DCOPY( N, B( 1, J ), 1, WORK( N+1 ), 1 )
CALL DGEMV( TRANS, N, N, -ONE, A, LDA, X( 1, J ), 1, ONE,
$ WORK( N+1 ), 1 )
*
* Compute componentwise relative backward error from formula
*
* max(i) ( abs(R(i)) / ( abs(op(A))*abs(X) + abs(B) )(i) )
*
* where abs(Z) is the componentwise absolute value of the matrix
* or vector Z. If the i-th component of the denominator is less
* than SAFE2, then SAFE1 is added to the i-th components of the
* numerator and denominator before dividing.
*
DO 30 I = 1, N
WORK( I ) = ABS( B( I, J ) )
30 CONTINUE
*
* Compute abs(op(A))*abs(X) + abs(B).
*
IF( NOTRAN ) THEN
DO 50 K = 1, N
XK = ABS( X( K, J ) )
DO 40 I = 1, N
WORK( I ) = WORK( I ) + ABS( A( I, K ) )*XK
40 CONTINUE
50 CONTINUE
ELSE
DO 70 K = 1, N
S = ZERO
DO 60 I = 1, N
S = S + ABS( A( I, K ) )*ABS( X( I, J ) )
60 CONTINUE
WORK( K ) = WORK( K ) + S
70 CONTINUE
END IF
S = ZERO
DO 80 I = 1, N
IF( WORK( I ).GT.SAFE2 ) THEN
S = MAX( S, ABS( WORK( N+I ) ) / WORK( I ) )
ELSE
S = MAX( S, ( ABS( WORK( N+I ) )+SAFE1 ) /
$ ( WORK( I )+SAFE1 ) )
END IF
80 CONTINUE
BERR( J ) = S
*
* Test stopping criterion. Continue iterating if
* 1) The residual BERR(J) is larger than machine epsilon, and
* 2) BERR(J) decreased by at least a factor of 2 during the
* last iteration, and
* 3) At most ITMAX iterations tried.
*
IF( BERR( J ).GT.EPS .AND. TWO*BERR( J ).LE.LSTRES .AND.
$ COUNT.LE.ITMAX ) THEN
*
* Update solution and try again.
*
CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK( N+1 ), N,
$ INFO )
CALL DAXPY( N, ONE, WORK( N+1 ), 1, X( 1, J ), 1 )
LSTRES = BERR( J )
COUNT = COUNT + 1
GO TO 20
END IF
*
* Bound error from formula
*
* norm(X - XTRUE) / norm(X) .le. FERR =
* norm( abs(inv(op(A)))*
* ( abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) ))) / norm(X)
*
* where
* norm(Z) is the magnitude of the largest component of Z
* inv(op(A)) is the inverse of op(A)
* abs(Z) is the componentwise absolute value of the matrix or
* vector Z
* NZ is the maximum number of nonzeros in any row of A, plus 1
* EPS is machine epsilon
*
* The i-th component of abs(R)+NZ*EPS*(abs(op(A))*abs(X)+abs(B))
* is incremented by SAFE1 if the i-th component of
* abs(op(A))*abs(X) + abs(B) is less than SAFE2.
*
* Use DLACON to estimate the infinity-norm of the matrix
* inv(op(A)) * diag(W),
* where W = abs(R) + NZ*EPS*( abs(op(A))*abs(X)+abs(B) )))
*
DO 90 I = 1, N
IF( WORK( I ).GT.SAFE2 ) THEN
WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I )
ELSE
WORK( I ) = ABS( WORK( N+I ) ) + NZ*EPS*WORK( I ) + SAFE1
END IF
90 CONTINUE
*
KASE = 0
100 CONTINUE
CALL DLACON( N, WORK( 2*N+1 ), WORK( N+1 ), IWORK, FERR( J ),
$ KASE )
IF( KASE.NE.0 ) THEN
IF( KASE.EQ.1 ) THEN
*
* Multiply by diag(W)*inv(op(A)**T).
*
CALL DGETRS( TRANST, N, 1, AF, LDAF, IPIV, WORK( N+1 ),
$ N, INFO )
DO 110 I = 1, N
WORK( N+I ) = WORK( I )*WORK( N+I )
110 CONTINUE
ELSE
*
* Multiply by inv(op(A))*diag(W).
*
DO 120 I = 1, N
WORK( N+I ) = WORK( I )*WORK( N+I )
120 CONTINUE
CALL DGETRS( TRANS, N, 1, AF, LDAF, IPIV, WORK( N+1 ), N,
$ INFO )
END IF
GO TO 100
END IF
*
* Normalize error.
*
LSTRES = ZERO
DO 130 I = 1, N
LSTRES = MAX( LSTRES, ABS( X( I, J ) ) )
130 CONTINUE
IF( LSTRES.NE.ZERO )
$ FERR( J ) = FERR( J ) / LSTRES
*
140 CONTINUE
*
RETURN
*
* End of DGERFS
*
END

148
ext/lapack/dtrtrs.f Normal file
View file

@ -0,0 +1,148 @@
SUBROUTINE DTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB,
$ INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* March 31, 1993
*
* .. Scalar Arguments ..
CHARACTER DIAG, TRANS, UPLO
INTEGER INFO, LDA, LDB, N, NRHS
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), B( LDB, * )
* ..
*
* Purpose
* =======
*
* DTRTRS solves a triangular system of the form
*
* A * X = B or A**T * X = B,
*
* where A is a triangular matrix of order N, and B is an N-by-NRHS
* matrix. A check is made to verify that A is nonsingular.
*
* Arguments
* =========
*
* UPLO (input) CHARACTER*1
* = 'U': A is upper triangular;
* = 'L': A is lower triangular.
*
* TRANS (input) CHARACTER*1
* Specifies the form of the system of equations:
* = 'N': A * X = B (No transpose)
* = 'T': A**T * X = B (Transpose)
* = 'C': A**H * X = B (Conjugate transpose = Transpose)
*
* DIAG (input) CHARACTER*1
* = 'N': A is non-unit triangular;
* = 'U': A is unit triangular.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* NRHS (input) INTEGER
* The number of right hand sides, i.e., the number of columns
* of the matrix B. NRHS >= 0.
*
* A (input) DOUBLE PRECISION array, dimension (LDA,N)
* The triangular matrix A. If UPLO = 'U', the leading N-by-N
* upper triangular part of the array A contains the upper
* triangular matrix, and the strictly lower triangular part of
* A is not referenced. If UPLO = 'L', the leading N-by-N lower
* triangular part of the array A contains the lower triangular
* matrix, and the strictly upper triangular part of A is not
* referenced. If DIAG = 'U', the diagonal elements of A are
* also not referenced and are assumed to be 1.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
* On entry, the right hand side matrix B.
* On exit, if INFO = 0, the solution matrix X.
*
* LDB (input) INTEGER
* The leading dimension of the array B. LDB >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value
* > 0: if INFO = i, the i-th diagonal element of A is zero,
* indicating that the matrix is singular and the solutions
* X have not been computed.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL NOUNIT
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DTRSM, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX
* ..
* .. Executable Statements ..
*
* Test the input parameters.
*
INFO = 0
NOUNIT = LSAME( DIAG, 'N' )
IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
INFO = -1
ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
$ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
INFO = -2
ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
INFO = -3
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( NRHS.LT.0 ) THEN
INFO = -5
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -7
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -9
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DTRTRS', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
*
* Check for singularity.
*
IF( NOUNIT ) THEN
DO 10 INFO = 1, N
IF( A( INFO, INFO ).EQ.ZERO )
$ RETURN
10 CONTINUE
END IF
INFO = 0
*
* Solve A * x = b or A' * x = b.
*
CALL DTRSM( 'Left', UPLO, TRANS, DIAG, N, NRHS, ONE, A, LDA, B,
$ LDB )
*
RETURN
*
* End of DTRTRS
*
END