From 144c6b23bde22d84756f75f5a367a3a85363b3ae Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Sat, 19 Feb 2011 00:08:29 +0000 Subject: [PATCH] Another delta on the dogleg implementation of the nonlinear solver. Added lapack routines for calculation of condition number and to fill out the QR factorization capability. --- Cantera/src/base/Array.h | 10 + Cantera/src/equil/vcs_solve.h | 4 +- Cantera/src/numerics/NonlinearSolver.cpp | 173 ++++++++-- Cantera/src/numerics/NonlinearSolver.h | 21 +- ext/f2c_lapack/Makefile.in | 4 + ext/f2c_lapack/dgecon.c | 228 ++++++++++++ ext/f2c_lapack/dgeequ.c | 297 ++++++++++++++++ ext/f2c_lapack/dgerfs.c | 420 +++++++++++++++++++++++ ext/f2c_lapack/dtrtrs.c | 187 ++++++++++ ext/lapack/Makefile.in | 4 + ext/lapack/dgecon.f | 181 ++++++++++ ext/lapack/dgeequ.f | 226 ++++++++++++ ext/lapack/dgerfs.f | 332 ++++++++++++++++++ ext/lapack/dtrtrs.f | 148 ++++++++ 14 files changed, 2200 insertions(+), 35 deletions(-) create mode 100644 ext/f2c_lapack/dgecon.c create mode 100644 ext/f2c_lapack/dgeequ.c create mode 100644 ext/f2c_lapack/dgerfs.c create mode 100644 ext/f2c_lapack/dtrtrs.c create mode 100644 ext/lapack/dgecon.f create mode 100644 ext/lapack/dgeequ.f create mode 100644 ext/lapack/dgerfs.f create mode 100644 ext/lapack/dtrtrs.f diff --git a/Cantera/src/base/Array.h b/Cantera/src/base/Array.h index a81864175..ae5e472df 100644 --- a/Cantera/src/base/Array.h +++ b/Cantera/src/base/Array.h @@ -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. diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index daf952d26..e0c52b12a 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -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 diff --git a/Cantera/src/numerics/NonlinearSolver.cpp b/Cantera/src/numerics/NonlinearSolver.cpp index 3f878ea6f..06a6b1bfe 100644 --- a/Cantera/src/numerics/NonlinearSolver.cpp +++ b/Cantera/src/numerics/NonlinearSolver.cpp @@ -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++) { + + + + + + + + + + + + + + } + + } //==================================================================================================================== /* * diff --git a/Cantera/src/numerics/NonlinearSolver.h b/Cantera/src/numerics/NonlinearSolver.h index ccda92128..47df6ab86 100644 --- a/Cantera/src/numerics/NonlinearSolver.h +++ b/Cantera/src/numerics/NonlinearSolver.h @@ -248,6 +248,7 @@ namespace Cantera { */ doublereal calcTrustDistance(std::vector 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 deltax_cp_; + std::vector deltaX_CP_; //! Newton Step - This is the newton step determined from the straight Jacobian /* @@ -822,7 +831,7 @@ namespace Cantera { std::vector Jd_; //! Vector of trust region values. - std::vector trustDeltaX_; + std::vector 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 *****************************************************************************************/ diff --git a/ext/f2c_lapack/Makefile.in b/ext/f2c_lapack/Makefile.in index f85c8c20a..68dff0ce8 100755 --- a/ext/f2c_lapack/Makefile.in +++ b/ext/f2c_lapack/Makefile.in @@ -100,6 +100,10 @@ dormqr.o \ drscl.o \ dtrtri.o \ dtrti2.o \ +dtrtrs.o \ +dgecon.o \ +dgeequ.o \ +dgerfs.o \ ieeeck.o \ ilaenv.o diff --git a/ext/f2c_lapack/dgecon.c b/ext/f2c_lapack/dgecon.c new file mode 100644 index 000000000..277e922af --- /dev/null +++ b/ext/f2c_lapack/dgecon.c @@ -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_ */ + diff --git a/ext/f2c_lapack/dgeequ.c b/ext/f2c_lapack/dgeequ.c new file mode 100644 index 000000000..125849ce7 --- /dev/null +++ b/ext/f2c_lapack/dgeequ.c @@ -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_ */ + diff --git a/ext/f2c_lapack/dgerfs.c b/ext/f2c_lapack/dgerfs.c new file mode 100644 index 000000000..5dcebe627 --- /dev/null +++ b/ext/f2c_lapack/dgerfs.c @@ -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_ */ + diff --git a/ext/f2c_lapack/dtrtrs.c b/ext/f2c_lapack/dtrtrs.c new file mode 100644 index 000000000..4c4035d8e --- /dev/null +++ b/ext/f2c_lapack/dtrtrs.c @@ -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_ */ + diff --git a/ext/lapack/Makefile.in b/ext/lapack/Makefile.in index 48b38590c..50a478cca 100755 --- a/ext/lapack/Makefile.in +++ b/ext/lapack/Makefile.in @@ -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) diff --git a/ext/lapack/dgecon.f b/ext/lapack/dgecon.f new file mode 100644 index 000000000..f6bd485f4 --- /dev/null +++ b/ext/lapack/dgecon.f @@ -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 diff --git a/ext/lapack/dgeequ.f b/ext/lapack/dgeequ.f new file mode 100644 index 000000000..a1be3c169 --- /dev/null +++ b/ext/lapack/dgeequ.f @@ -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 diff --git a/ext/lapack/dgerfs.f b/ext/lapack/dgerfs.f new file mode 100644 index 000000000..aa6c13415 --- /dev/null +++ b/ext/lapack/dgerfs.f @@ -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 diff --git a/ext/lapack/dtrtrs.f b/ext/lapack/dtrtrs.f new file mode 100644 index 000000000..c1b4c5c4c --- /dev/null +++ b/ext/lapack/dtrtrs.f @@ -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