diff --git a/lapack-netlib/SRC/cggglm.c b/lapack-netlib/SRC/cggglm.c index e04562eb1c..b72d182c4c 100644 --- a/lapack-netlib/SRC/cggglm.c +++ b/lapack-netlib/SRC/cggglm.c @@ -216,7 +216,11 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +#define myexp_(w) my_expfunc(w) +static int my_expfunc(float *x) {int e; (void)frexp(*x,&e); return e;} +static float myscale_(float *w, integer *i) {float myscale=*w*pow(2.,*(double*)i); return(myscale);} +static float myhuge_(float *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -226,13 +230,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static complex c_b2 = {1.f,0.f}; static integer c__1 = 1; static integer c_n1 = -1; +static real c_b26 = 0.f; +static integer c__0 = 0; +static real c_b30 = 1.f; /* > \brief \b CGGGLM */ @@ -241,7 +246,6 @@ static integer c_n1 = -1; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download CGGGLM + dependencies */ /* > */ @@ -252,7 +256,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -296,6 +299,16 @@ f"> */ /* > x */ /* > */ /* > where inv(B) denotes the inverse of B. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The CTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized QR */ +/* > factorization of the pair (A, B) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -394,12 +407,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with A in the */ -/* > generalized QR factorization of the pair (A, B) is */ +/* > generalized QR factorization of the pair (A, B) is exactly */ /* > singular, so that rank(A) < M; the least squares */ /* > solution could not be computed. */ /* > = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal */ /* > factor T associated with B in the generalized QR */ -/* > factorization of the pair (A, B) is singular, so that */ +/* > factorization of the pair (A, B) is exactly singular, so that */ /* > rank( A B ) < N; the least squares solution could not */ /* > be computed. */ /* > \endverbatim */ @@ -412,9 +425,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complexOTHEReigen */ +/* > \ingroup ggglm */ /* ===================================================================== */ /* Subroutine */ void cggglm_(integer *n, integer *m, integer *p, complex *a, @@ -423,37 +434,49 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; - complex q__1={0.,0.}; + real r__1, r__2; + complex q__1; /* Local variables */ - integer lopt, i__; + integer i__, ia, ib, ka, id, kb, kd, nb, np, nb1, nb2, nb3, nb4, ibig; + real anrm, bnrm, dnrm; + integer isml, lopt; extern /* Subroutine */ void cgemv_(char *, integer *, integer *, complex * , complex *, integer *, complex *, integer *, complex *, complex * , integer *), ccopy_(integer *, complex *, integer *, complex *, integer *); - integer nb, np; - extern /* Subroutine */ void cggqrf_(integer *, integer *, integer *, - complex *, integer *, complex *, complex *, integer *, complex *, - complex *, integer *, integer *); - extern void xerbla_(char *, integer *, ftnlen); + real rwork[1]; + extern real clange_(char *, integer *, integer *, complex *, integer *, + real *); + extern /* Subroutine */ void clascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, complex *, integer *, integer *); + extern real slamch_(char *); + extern /* Subroutine */ void csscal_(integer *, real *, complex *, integer + *), cggqrf_(integer *, integer *, integer *, complex *, integer *, + complex *, complex *, integer *, complex *, complex *, integer *, + integer *), xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - integer lwkmin, nb1, nb2, nb3, nb4; + real bignum; + //extern real myhuge_(real *); + integer lwkmin; extern /* Subroutine */ void cunmqr_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *), cunmrq_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *); + real smlnum; integer lwkopt; logical lquery; extern /* Subroutine */ void ctrtrs_(char *, char *, char *, integer *, integer *, complex *, integer *, complex *, integer *, integer *); + extern real sroundup_lwork__(integer *); + //, myscale_(real *, integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* =================================================================== */ @@ -463,10 +486,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --d__; --x; @@ -510,7 +533,8 @@ f"> */ lwkmin = *m + *n + *p; lwkopt = *m + np + f2cmax(*n,*p) * nb; } - work[1].r = (real) lwkopt, work[1].i = 0.f; + r__1 = sroundup_lwork__(&lwkopt); + work[1].r = r__1, work[1].i = 0.f; if (*lwork < lwkmin && ! lquery) { *info = -12; @@ -541,6 +565,78 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = slamch_("S") / slamch_("P"); + bignum = 1.f / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B and d independently by powers of two 2**KA, 2**KB and */ +/* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the */ +/* values whose MYEXP lies in [ISML,IBIG]. x is then scaled by */ +/* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is */ +/* exact. A norm that is zero, infinite or NaN takes no part. */ +/* CLASCL is called with both endpoints at or above one so that */ +/* forming the factor raises no underflow. */ + + anrm = clange_("M", n, m, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0.f && anrm <= myhuge_(&c_b26)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b30, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b30, &i__3); + clascl_("G", &c__0, &c__0, &r__1, &r__2, n, m, &a[a_offset], lda, + info); + } + + bnrm = clange_("M", n, p, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0.f && bnrm <= myhuge_(&c_b26)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b30, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b30, &i__3); + clascl_("G", &c__0, &c__0, &r__1, &r__2, n, p, &b[b_offset], ldb, + info); + } + + dnrm = clange_("M", n, &c__1, &d__[1], n, rwork); + kd = 0; + if (dnrm > 0.f && dnrm <= myhuge_(&c_b26)) { + id = myexp_(&dnrm); + if (id < isml) { + kd = isml - id; + } else if (id > ibig) { + kd = ibig - id; + } + } + if (kd != 0) { + r__1 = myscale_(&c_b30, &kd); + csscal_(n, &r__1, &d__[1], &c__1); + } + /* Compute the GQR factorization of matrices A and B: */ /* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M */ @@ -554,7 +650,7 @@ f"> */ cggqrf_(n, m, p, &a[a_offset], lda, &work[1], &b[b_offset], ldb, &work[*m + 1], &work[*m + np + 1], &i__1, info); i__1 = *m + np + 1; - lopt = work[i__1].r; + lopt = (integer) work[i__1].r; /* Update left-hand-side vector d = Q**H*d = ( d1 ) M */ /* ( d2 ) N-M */ @@ -598,7 +694,7 @@ f"> */ /* Update d1 = d1 - T12*y2 */ i__1 = *n - *m; - q__1.r = -1.f, q__1.i = 0.f; + q__1.r = -1.f, q__1.i = -0.f; cgemv_("No transpose", m, &i__1, &q__1, &b[(*m + *p - *n + 1) * b_dim1 + 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b2, &d__[1], &c__1); @@ -627,11 +723,26 @@ f"> */ cunmrq_("Left", "Conjugate transpose", p, &c__1, &np, &b[f2cmax(i__1,i__2) + b_dim1], ldb, &work[*m + 1], &y[1], &i__3, &work[*m + np + 1], & i__4, info); + +/* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) */ + + if (ka != kd) { + i__1 = ka - kd; + r__1 = myscale_(&c_b30, &i__1); + csscal_(m, &r__1, &x[1], &c__1); + } + if (kb != kd) { + i__1 = kb - kd; + r__1 = myscale_(&c_b30, &i__1); + csscal_(p, &r__1, &y[1], &c__1); + } + /* Computing MAX */ i__4 = *m + np + 1; i__2 = lopt, i__3 = (integer) work[i__4].r; i__1 = *m + np + f2cmax(i__2,i__3); - work[1].r = (real) i__1, work[1].i = 0.f; + q__1.r = (real) i__1, q__1.i = 0.f; + work[1].r = q__1.r, work[1].i = q__1.i; return; diff --git a/lapack-netlib/SRC/cgglse.c b/lapack-netlib/SRC/cgglse.c index f62f7e44af..3b06184221 100644 --- a/lapack-netlib/SRC/cgglse.c +++ b/lapack-netlib/SRC/cgglse.c @@ -216,7 +216,13 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +#define myexp_(w) my_expfunc(w) +static int my_expfunc(float *x) {int e; (void)frexp(*x,&e); return e;} +static float myscale_(float *w, integer *i) {return (*w*pow(2.,*(float*)i));} + +static float myhuge_(float *w) {return(HUGE_VAL);} +static int myihuge_(int *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -226,13 +232,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static complex c_b1 = {1.f,0.f}; static integer c__1 = 1; static integer c_n1 = -1; +static real c_b25 = 0.f; +static integer c__0 = 0; +static real c_b29 = 1.f; /* > \brief CGGLSE solves overdetermined or underdetermined systems for OTHER matrices */ @@ -241,7 +248,6 @@ static integer c_n1 = -1; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download CGGLSE + dependencies */ /* > */ @@ -252,7 +258,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -287,6 +292,16 @@ f"> */ /* > matrices (B, A) given by */ /* > */ /* > B = (0 R)*Q, A = Z*T*Q. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The CTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized RQ */ +/* > factorization of the pair (B, A) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -388,12 +403,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with B in the */ -/* > generalized RQ factorization of the pair (B, A) is */ +/* > generalized RQ factorization of the pair (B, A) is exactly */ /* > singular, so that rank(B) < P; the least squares */ /* > solution could not be computed. */ /* > = 2: the (N-P) by (N-P) part of the upper trapezoidal factor */ /* > T associated with A in the generalized RQ factorization */ -/* > of the pair (B, A) is singular, so that */ +/* > of the pair (B, A) is exactly singular, so that */ /* > rank( (A) ) < N; the least squares solution could not */ /* > ( (B) ) */ /* > be computed. */ @@ -407,9 +422,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complexOTHERsolve */ +/* > \ingroup gglse */ /* ===================================================================== */ /* Subroutine */ void cgglse_(integer *m, integer *n, integer *p, complex *a, @@ -418,39 +431,51 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; - complex q__1={0.,0.}; + real r__1, r__2; + complex q__1; /* Local variables */ - integer lopt; + integer ia, ib, ic, id, ka, kb, nb, mn, iu, kt, nr, nb1, nb2, nb3, nb4, + ibig; + real anrm, bnrm, cnrm, dnrm; + integer isml, lopt; extern /* Subroutine */ void cgemv_(char *, integer *, integer *, complex * , complex *, integer *, complex *, integer *, complex *, complex * , integer *), ccopy_(integer *, complex *, integer *, complex *, integer *), caxpy_(integer *, complex *, complex *, integer *, complex *, integer *), ctrmv_(char *, char *, char *, integer *, complex *, integer *, complex *, integer *); - integer nb, mn, nr; - extern /* Subroutine */ void cggrqf_(integer *, integer *, integer *, - complex *, integer *, complex *, complex *, integer *, complex *, - complex *, integer *, integer *); - extern void xerbla_(char *, integer *, ftnlen); + extern integer myexp_(real *); + real rwork[1]; + extern real clange_(char *, integer *, integer *, complex *, integer *, + real *); + extern /* Subroutine */ void clascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, complex *, integer *, integer *); + extern real slamch_(char *); + extern /* Subroutine */ void csscal_(integer *, real *, complex *, integer + *), cggrqf_(integer *, integer *, integer *, complex *, integer *, + complex *, complex *, integer *, complex *, complex *, integer *, + integer *), xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - integer lwkmin, nb1, nb2, nb3, nb4; + real bignum; + integer lwkmin; extern /* Subroutine */ void cunmqr_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *), cunmrq_(char *, char *, integer *, integer *, integer *, complex *, integer *, complex *, complex *, integer *, complex *, integer *, integer *); + real smlnum; integer lwkopt; logical lquery; extern /* Subroutine */ void ctrtrs_(char *, char *, char *, integer *, integer *, complex *, integer *, complex *, integer *, integer *); + extern real sroundup_lwork__(integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* ===================================================================== */ @@ -460,10 +485,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --c__; --d__; @@ -507,7 +532,8 @@ f"> */ lwkmin = *m + *n + *p; lwkopt = *p + mn + f2cmax(*m,*n) * nb; } - work[1].r = (real) lwkopt, work[1].i = 0.f; + r__1 = sroundup_lwork__(&lwkopt); + work[1].r = r__1, work[1].i = 0.f; if (*lwork < lwkmin && ! lquery) { *info = -12; @@ -528,6 +554,96 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = slamch_("S") / slamch_("P"); + bignum = 1.f / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B, c and d by powers of two so that their largest */ +/* entries lie in [SMLNUM,BIGNUM), the values whose MYEXP lies */ +/* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and */ +/* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which */ +/* keeps the constraint and the residual consistent and scales x by */ +/* 2**KT; KT brings the larger of the two right-hand sides into */ +/* range. Scaling by a power of two is exact. A norm that is zero, */ +/* infinite or NaN takes no part. CLASCL is called with both */ +/* endpoints at or above one so that forming the factor raises no */ +/* underflow. */ + + anrm = clange_("M", m, n, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0.f && anrm <= myhuge_(&c_b25)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b29, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b29, &i__3); + clascl_("G", &c__0, &c__0, &r__1, &r__2, m, n, &a[a_offset], lda, + info); + } + + bnrm = clange_("M", p, n, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0.f && bnrm <= myhuge_(&c_b25)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b29, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b29, &i__3); + clascl_("G", &c__0, &c__0, &r__1, &r__2, p, n, &b[b_offset], ldb, + info); + } + + i__1 = f2cmax(1,*m); + cnrm = clange_("M", m, &c__1, &c__[1], &i__1, rwork); + i__1 = f2cmax(1,*p); + dnrm = clange_("M", p, &c__1, &d__[1], &i__1, rwork); + ic = -myihuge_(&c__0); + if (cnrm > 0.f && cnrm <= myhuge_(&c_b25)) { + ic = myexp_(&cnrm) + ka; + } + id = -myihuge_(&c__0); + if (dnrm > 0.f && dnrm <= myhuge_(&c_b25)) { + id = myexp_(&dnrm) + kb; + } + iu = f2cmax(ic,id); + kt = 0; + if (iu > ibig) { + kt = ibig - iu; + } else if (iu < isml && (real) iu != -myihuge_(&c__0)) { + kt = isml - iu; + } + if (ka + kt != 0) { + i__1 = ka + kt; + r__1 = myscale_(&c_b29, &i__1); + csscal_(m, &r__1, &c__[1], &c__1); + } + if (kb + kt != 0) { + i__1 = kb + kt; + r__1 = myscale_(&c_b29, &i__1); + csscal_(p, &r__1, &d__[1], &c__1); + } + /* Compute the GRQ factorization of matrices B and A: */ /* B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P */ @@ -541,7 +657,7 @@ f"> */ cggrqf_(p, m, n, &b[b_offset], ldb, &work[1], &a[a_offset], lda, &work[*p + 1], &work[*p + mn + 1], &i__1, info); i__1 = *p + mn + 1; - lopt = work[i__1].r; + lopt = (integer) work[i__1].r; /* Update c = Z**H *c = ( c1 ) N-P */ /* ( c2 ) M+P-N */ @@ -573,7 +689,7 @@ f"> */ /* Update c1 */ i__1 = *n - *p; - q__1.r = -1.f, q__1.i = 0.f; + q__1.r = -1.f, q__1.i = -0.f; cgemv_("No transpose", &i__1, p, &q__1, &a[(*n - *p + 1) * a_dim1 + 1] , lda, &d__[1], &c__1, &c_b1, &c__[1], &c__1); } @@ -603,7 +719,7 @@ f"> */ nr = *m + *p - *n; if (nr > 0) { i__1 = *n - *m; - q__1.r = -1.f, q__1.i = 0.f; + q__1.r = -1.f, q__1.i = -0.f; cgemv_("No transpose", &nr, &i__1, &q__1, &a[*n - *p + 1 + (*m + 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b1, &c__[*n - * p + 1], &c__1); @@ -614,7 +730,7 @@ f"> */ if (nr > 0) { ctrmv_("Upper", "No transpose", "Non unit", &nr, &a[*n - *p + 1 + (*n - *p + 1) * a_dim1], lda, &d__[1], &c__1); - q__1.r = -1.f, q__1.i = 0.f; + q__1.r = -1.f, q__1.i = -0.f; caxpy_(&nr, &q__1, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); } @@ -623,11 +739,28 @@ f"> */ i__1 = *lwork - *p - mn; cunmrq_("Left", "Conjugate Transpose", n, &c__1, p, &b[b_offset], ldb, & work[1], &x[1], n, &work[*p + mn + 1], &i__1, info); + +/* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) */ +/* carries 2**(KA+KT) */ + + if (kt != 0) { + i__1 = -kt; + r__1 = myscale_(&c_b29, &i__1); + csscal_(n, &r__1, &x[1], &c__1); + } + if (ka + kt != 0 && *m + *p > *n) { + i__1 = *m + *p - *n; + i__2 = -ka - kt; + r__1 = myscale_(&c_b29, &i__2); + csscal_(&i__1, &r__1, &c__[*n - *p + 1], &c__1); + } + /* Computing MAX */ i__4 = *p + mn + 1; i__2 = lopt, i__3 = (integer) work[i__4].r; i__1 = *p + mn + f2cmax(i__2,i__3); - work[1].r = (real) i__1, work[1].i = 0.f; + q__1.r = (real) i__1, q__1.i = 0.f; + work[1].r = q__1.r, work[1].i = q__1.i; return; diff --git a/lapack-netlib/SRC/cgglse.f b/lapack-netlib/SRC/cgglse.f index b1c5623858..1a24122bb8 100644 --- a/lapack-netlib/SRC/cgglse.f +++ b/lapack-netlib/SRC/cgglse.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download CGGLSE + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -52,6 +50,16 @@ *> matrices (B, A) given by *> *> B = (0 R)*Q, A = Z*T*Q. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The CTRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized RQ +*> factorization of the pair (B, A) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -153,12 +161,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with B in the -*> generalized RQ factorization of the pair (B, A) is +*> generalized RQ factorization of the pair (B, A) is exactly *> singular, so that rank(B) < P; the least squares *> solution could not be computed. *> = 2: the (N-P) by (N-P) part of the upper trapezoidal factor *> T associated with A in the generalized RQ factorization -*> of the pair (B, A) is singular, so that +*> of the pair (B, A) is exactly singular, so that *> rank( (A) ) < N; the least squares solution could not *> ( (B) ) *> be computed. @@ -175,8 +183,10 @@ *> \ingroup gglse * * ===================================================================== - SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, + SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -193,6 +203,8 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * ===================================================================== * * .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) COMPLEX CONE PARAMETER ( CONE = ( 1.0E+0, 0.0E+0 ) ) * .. @@ -200,18 +212,24 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, LOGICAL LQUERY INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, $ NB4, NR + INTEGER IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT + REAL ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + REAL RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL CAXPY, CCOPY, CGEMV, CGGRQF, CTRMV, CTRTRS, - $ CUNMQR, CUNMRQ, XERBLA + EXTERNAL CAXPY, CCOPY, CGEMV, CGGRQF, CLASCL, CSSCAL, + $ CTRMV, CTRTRS, CUNMQR, CUNMRQ, XERBLA * .. * .. External Functions .. INTEGER ILAENV REAL SROUNDUP_LWORK - EXTERNAL ILAENV, SROUNDUP_LWORK + REAL SLAMCH, CLANGE + EXTERNAL CLANGE, ILAENV, SLAMCH, SROUNDUP_LWORK * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -266,6 +284,72 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( N.EQ.0 ) $ RETURN * +* Get machine parameters +* + SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B, c and d by powers of two so that their largest +* entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies +* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and +* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which +* keeps the constraint and the residual consistent and scales x by +* 2**KT; KT brings the larger of the two right-hand sides into +* range. Scaling by a power of two is exact. A norm that is zero, +* infinite or NaN takes no part. CLASCL is called with both +* endpoints at or above one so that forming the factor raises no +* underflow. +* + ANRM = CLANGE( 'M', M, N, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO ) +* + BNRM = CLANGE( 'M', P, N, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL CLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO ) +* + CNRM = CLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK ) + DNRM = CLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK ) + IC = -HUGE( 0 ) + IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) ) + $ IC = EXPONENT( CNRM ) + KA + ID = -HUGE( 0 ) + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) + $ ID = EXPONENT( DNRM ) + KB + IU = MAX( IC, ID ) + KT = 0 + IF( IU.GT.IBIG ) THEN + KT = IBIG - IU + ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN + KT = ISML - IU + END IF + IF( KA+KT.NE.0 ) + $ CALL CSSCAL( M, SCALE( ONE, KA+KT ), C, 1 ) + IF( KB+KT.NE.0 ) + $ CALL CSSCAL( P, SCALE( ONE, KB+KT ), D, 1 ) +* * Compute the GRQ factorization of matrices B and A: * * B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P @@ -304,7 +388,8 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Update c1 * - CALL CGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA, + CALL CGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), + $ LDA, $ D, 1, CONE, C, 1 ) END IF * @@ -329,7 +414,8 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( M.LT.N ) THEN NR = M + P - N IF( NR.GT.0 ) - $ CALL CGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ), + $ CALL CGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, + $ M+1 ), $ LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 ) ELSE NR = P @@ -344,7 +430,16 @@ SUBROUTINE CGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * CALL CUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB, $ WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) - WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) +* +* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) +* carries 2**(KA+KT) +* + IF( KT.NE.0 ) + $ CALL CSSCAL( N, SCALE( ONE, -KT ), X, 1 ) + IF( KA+KT.NE.0 .AND. M+P.GT.N ) + $ CALL CSSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 ) +* + WORK( 1 ) = CMPLX( P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) ) * RETURN * diff --git a/lapack-netlib/SRC/dggglm.c b/lapack-netlib/SRC/dggglm.c index 4c83b89da4..92642c926b 100644 --- a/lapack-netlib/SRC/dggglm.c +++ b/lapack-netlib/SRC/dggglm.c @@ -216,7 +216,11 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(double *x) {int e; (void)frexp(*x,&e); return e;} +static double myscale_(double *w, integer *i) {double myscale=*w*pow(2.,*(double*)i); return(myscale);} +static double myhuge_(double *w) {return(HUGE_VAL);} + + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -226,14 +230,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c_n1 = -1; -static doublereal c_b32 = -1.; -static doublereal c_b34 = 1.; +static doublereal c_b24 = 0.; +static integer c__0 = 0; +static doublereal c_b28 = 1.; +static doublereal c_b53 = -1.; /* > \brief \b DGGGLM */ @@ -242,7 +246,6 @@ static doublereal c_b34 = 1.; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download DGGGLM + dependencies */ /* > */ @@ -253,7 +256,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -297,6 +299,16 @@ f"> */ /* > x */ /* > */ /* > where inv(B) denotes the inverse of B. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The DTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized QR */ +/* > factorization of the pair (A, B) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -395,12 +407,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with A in the */ -/* > generalized QR factorization of the pair (A, B) is */ +/* > generalized QR factorization of the pair (A, B) is exactly */ /* > singular, so that rank(A) < M; the least squares */ /* > solution could not be computed. */ /* > = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal */ /* > factor T associated with B in the generalized QR */ -/* > factorization of the pair (A, B) is singular, so that */ +/* > factorization of the pair (A, B) is exactly singular, so that */ /* > rank( A B ) < N; the least squares solution could not */ /* > be computed. */ /* > \endverbatim */ @@ -413,9 +425,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup doubleOTHEReigen */ +/* > \ingroup ggglm */ /* ===================================================================== */ /* Subroutine */ void dggglm_(integer *n, integer *m, integer *p, doublereal * @@ -425,27 +435,36 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; + real r__1, r__2; /* Local variables */ - integer lopt, i__; - extern /* Subroutine */ void dgemv_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, doublereal *, integer *, - doublereal *, doublereal *, integer *), dcopy_(integer *, - doublereal *, integer *, doublereal *, integer *); - integer nb, np; - extern /* Subroutine */ void dggqrf_(integer *, integer *, integer *, - doublereal *, integer *, doublereal *, doublereal *, integer *, - doublereal *, doublereal *, integer *, integer *); - extern void xerbla_(char *, integer *, ftnlen); + integer i__, ia, ib, ka, id, kb, kd, nb, np, nb1, nb2, nb3, nb4, ibig; + doublereal anrm, bnrm, dnrm; + integer isml, lopt; + extern /* Subroutine */ void dscal_(integer *, real *, doublereal *, + integer *), dgemv_(char *, integer *, integer *, doublereal *, + doublereal *, integer *, doublereal *, integer *, doublereal *, + doublereal *, integer *), dcopy_(integer *, doublereal *, + integer *, doublereal *, integer *); + doublereal rwork[1]; + extern doublereal dlamch_(char *), dlange_(char *, integer *, + integer *, doublereal *, integer *, doublereal *); + extern /* Subroutine */ void dlascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, doublereal *, integer *, integer *), + dggqrf_(integer *, integer *, integer *, doublereal *, + integer *, doublereal *, doublereal *, integer *, doublereal *, + doublereal *, integer *, integer *), xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - integer lwkmin, nb1, nb2, nb3, nb4; + doublereal bignum; + integer lwkmin; extern /* Subroutine */ void dormqr_(char *, char *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *), dormrq_(char *, char *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *); + doublereal smlnum; integer lwkopt; logical lquery; extern /* Subroutine */ void dtrtrs_(char *, char *, char *, integer *, @@ -453,10 +472,9 @@ f"> */ integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* =================================================================== */ @@ -466,10 +484,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --d__; --x; @@ -542,6 +560,78 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = dlamch_("S") / dlamch_("P"); + bignum = 1. / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B and d independently by powers of two 2**KA, 2**KB and */ +/* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the */ +/* values whose MYEXP lies in [ISML,IBIG]. x is then scaled by */ +/* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is */ +/* exact. A norm that is zero, infinite or NaN takes no part. */ +/* DLASCL is called with both endpoints at or above one so that */ +/* forming the factor raises no underflow. */ + + anrm = dlange_("M", n, m, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0. && anrm <= myhuge_(&c_b24)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b28, &i__3); + dlascl_("G", &c__0, &c__0, &r__1, &r__2, n, m, &a[a_offset], lda, + info); + } + + bnrm = dlange_("M", n, p, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0. && bnrm <= myhuge_(&c_b24)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b28, &i__3); + dlascl_("G", &c__0, &c__0, &r__1, &r__2, n, p, &b[b_offset], ldb, + info); + } + + dnrm = dlange_("M", n, &c__1, &d__[1], n, rwork); + kd = 0; + if (dnrm > 0. && dnrm <= myhuge_(&c_b24)) { + id = myexp_(&dnrm); + if (id < isml) { + kd = isml - id; + } else if (id > ibig) { + kd = ibig - id; + } + } + if (kd != 0) { + r__1 = myscale_(&c_b28, &kd); + dscal_(n, &r__1, &d__[1], &c__1); + } + /* Compute the GQR factorization of matrices A and B: */ /* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M */ @@ -596,8 +686,8 @@ f"> */ /* Update d1 = d1 - T12*y2 */ i__1 = *n - *m; - dgemv_("No transpose", m, &i__1, &c_b32, &b[(*m + *p - *n + 1) * b_dim1 + - 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b34, &d__[1], &c__1); + dgemv_("No transpose", m, &i__1, &c_b53, &b[(*m + *p - *n + 1) * b_dim1 + + 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b28, &d__[1], &c__1); /* Solve triangular system: R11*x = d1 */ @@ -623,6 +713,20 @@ f"> */ i__4 = *lwork - *m - np; dormrq_("Left", "Transpose", p, &c__1, &np, &b[f2cmax(i__1,i__2) + b_dim1], ldb, &work[*m + 1], &y[1], &i__3, &work[*m + np + 1], &i__4, info); + +/* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) */ + + if (ka != kd) { + i__1 = ka - kd; + r__1 = myscale_(&c_b28, &i__1); + dscal_(m, &r__1, &x[1], &c__1); + } + if (kb != kd) { + i__1 = kb - kd; + r__1 = myscale_(&c_b28, &i__1); + dscal_(p, &r__1, &y[1], &c__1); + } + /* Computing MAX */ i__1 = lopt, i__2 = (integer) work[*m + np + 1]; work[1] = (doublereal) (*m + np + f2cmax(i__1,i__2)); diff --git a/lapack-netlib/SRC/dggglm.f b/lapack-netlib/SRC/dggglm.f index ae0f0e908c..0a46b6f618 100644 --- a/lapack-netlib/SRC/dggglm.f +++ b/lapack-netlib/SRC/dggglm.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DGGGLM + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -61,6 +59,16 @@ *> x *> *> where inv(B) denotes the inverse of B. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The DTRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized QR +*> factorization of the pair (A, B) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -159,12 +167,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with A in the -*> generalized QR factorization of the pair (A, B) is +*> generalized QR factorization of the pair (A, B) is exactly *> singular, so that rank(A) < M; the least squares *> solution could not be computed. *> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal *> factor T associated with B in the generalized QR -*> factorization of the pair (A, B) is singular, so that +*> factorization of the pair (A, B) is exactly singular, so that *> rank( A B ) < N; the least squares solution could not *> be computed. *> \endverbatim @@ -177,11 +185,13 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHEReigen +*> \ingroup ggglm * * ===================================================================== - SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, + SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -205,17 +215,23 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, LOGICAL LQUERY INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, $ NB4, NP + INTEGER IA, IB, IBIG, ID, ISML, KA, KB, KD + DOUBLE PRECISION ANRM, BIGNUM, BNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL DCOPY, DGEMV, DGGQRF, DORMQR, DORMRQ, DTRTRS, - $ XERBLA + EXTERNAL DCOPY, DGEMV, DGGQRF, DLASCL, DORMQR, DORMRQ, + $ DSCAL, DTRTRS, XERBLA * .. * .. External Functions .. INTEGER ILAENV - EXTERNAL ILAENV + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL DLAMCH, DLANGE, ILAENV * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -277,6 +293,62 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, RETURN END IF * +* Get machine parameters +* + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B and d independently by powers of two 2**KA, 2**KB and +* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the +* values whose EXPONENT lies in [ISML,IBIG]. x is then scaled by +* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is +* exact. A norm that is zero, infinite or NaN takes no part. +* DLASCL is called with both endpoints at or above one so that +* forming the factor raises no underflow. +* + ANRM = DLANGE( 'M', N, M, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO ) +* + BNRM = DLANGE( 'M', N, P, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO ) +* + DNRM = DLANGE( 'M', N, 1, D, N, RWORK ) + KD = 0 + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN + ID = EXPONENT( DNRM ) + IF( ID.LT.ISML ) THEN + KD = ISML - ID + ELSE IF( ID.GT.IBIG ) THEN + KD = IBIG - ID + END IF + END IF + IF( KD.NE.0 ) + $ CALL DSCAL( N, SCALE( ONE, KD ), D, 1 ) +* * Compute the GQR factorization of matrices A and B: * * Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M @@ -325,7 +397,8 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, * Solve triangular system: R11*x = d1 * IF( M.GT.0 ) THEN - CALL DTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, + CALL DTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, + $ LDA, $ D, M, INFO ) * IF( INFO.GT.0 ) THEN @@ -343,6 +416,14 @@ SUBROUTINE DGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, CALL DORMRQ( 'Left', 'Transpose', P, 1, NP, $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) +* +* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) +* + IF( KA.NE.KD ) + $ CALL DSCAL( M, SCALE( ONE, KA-KD ), X, 1 ) + IF( KB.NE.KD ) + $ CALL DSCAL( P, SCALE( ONE, KB-KD ), Y, 1 ) +* WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) * RETURN diff --git a/lapack-netlib/SRC/dgglse.c b/lapack-netlib/SRC/dgglse.c index 213de2e588..d73f257200 100644 --- a/lapack-netlib/SRC/dgglse.c +++ b/lapack-netlib/SRC/dgglse.c @@ -216,7 +216,10 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(double *x) {int e; (void)frexp(*x,&e); return e;} +static double myscale_(double *w, integer *i) {double myscale=*w*pow(2.,*(double*)i); return(myscale);} +static double myhuge_(double *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -226,14 +229,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c_n1 = -1; -static doublereal c_b31 = -1.; -static doublereal c_b33 = 1.; +static doublereal c_b24 = 0.; +static integer c__0 = 0; +static doublereal c_b28 = 1.; +static doublereal c_b60 = -1.; /* > \brief DGGLSE solves overdetermined or underdetermined systems for OTHER matrices */ @@ -242,7 +245,6 @@ static doublereal c_b33 = 1.; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download DGGLSE + dependencies */ /* > */ @@ -253,7 +255,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -288,6 +289,16 @@ f"> */ /* > matrices (B, A) given by */ /* > */ /* > B = (0 R)*Q, A = Z*T*Q. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The DTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized RQ */ +/* > factorization of the pair (B, A) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -389,12 +400,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with B in the */ -/* > generalized RQ factorization of the pair (B, A) is */ +/* > generalized RQ factorization of the pair (B, A) is exactly */ /* > singular, so that rank(B) < P; the least squares */ /* > solution could not be computed. */ /* > = 2: the (N-P) by (N-P) part of the upper trapezoidal factor */ /* > T associated with A in the generalized RQ factorization */ -/* > of the pair (B, A) is singular, so that */ +/* > of the pair (B, A) is exactly singular, so that */ /* > rank( (A) ) < N; the least squares solution could not */ /* > ( (B) ) */ /* > be computed. */ @@ -408,9 +419,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup doubleOTHERsolve */ +/* > \ingroup gglse */ /* ===================================================================== */ /* Subroutine */ void dgglse_(integer *m, integer *n, integer *p, doublereal * @@ -419,31 +428,41 @@ f"> */ integer *info) { /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; + integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3; + real r__1, r__2; /* Local variables */ - integer lopt; - extern /* Subroutine */ void dgemv_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, doublereal *, integer *, - doublereal *, doublereal *, integer *), dcopy_(integer *, - doublereal *, integer *, doublereal *, integer *), daxpy_(integer - *, doublereal *, doublereal *, integer *, doublereal *, integer *) - , dtrmv_(char *, char *, char *, integer *, doublereal *, integer - *, doublereal *, integer *); - integer nb, mn, nr; - extern /* Subroutine */ void dggrqf_(integer *, integer *, integer *, - doublereal *, integer *, doublereal *, doublereal *, integer *, - doublereal *, doublereal *, integer *, integer *); - extern void xerbla_(char *, integer *, ftnlen); + integer ia, ib, ic, id, ka, kb, nb, mn, iu, kt, nr, nb1, nb2, nb3, nb4, + ibig; + doublereal anrm, bnrm, cnrm, dnrm; + integer isml, lopt; + extern /* Subroutine */ void dscal_(integer *, real *, doublereal *, + integer *), dgemv_(char *, integer *, integer *, doublereal *, + doublereal *, integer *, doublereal *, integer *, doublereal *, + doublereal *, integer *), dcopy_(integer *, doublereal *, + integer *, doublereal *, integer *), daxpy_(integer *, doublereal + *, doublereal *, integer *, doublereal *, integer *), dtrmv_(char + *, char *, char *, integer *, doublereal *, integer *, doublereal + *, integer *); + doublereal rwork[1]; + extern doublereal dlamch_(char *), dlange_(char *, integer *, + integer *, doublereal *, integer *, doublereal *); + extern /* Subroutine */ void dlascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, doublereal *, integer *, integer *), + dggrqf_(integer *, integer *, integer *, doublereal *, + integer *, doublereal *, doublereal *, integer *, doublereal *, + doublereal *, integer *, integer *), xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - integer lwkmin, nb1, nb2, nb3, nb4; + doublereal bignum; + integer lwkmin; extern /* Subroutine */ void dormqr_(char *, char *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *), dormrq_(char *, char *, integer *, integer *, integer *, doublereal *, integer *, doublereal *, doublereal *, integer *, doublereal *, integer *, integer *); + doublereal smlnum; integer lwkopt; logical lquery; extern /* Subroutine */ void dtrtrs_(char *, char *, char *, integer *, @@ -451,10 +470,9 @@ f"> */ integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* ===================================================================== */ @@ -464,10 +482,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --c__; --d__; @@ -532,6 +550,96 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = dlamch_("S") / dlamch_("P"); + bignum = 1. / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B, c and d by powers of two so that their largest */ +/* entries lie in [SMLNUM,BIGNUM), the values whose MYEXP lies */ +/* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and */ +/* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which */ +/* keeps the constraint and the residual consistent and scales x by */ +/* 2**KT; KT brings the larger of the two right-hand sides into */ +/* range. Scaling by a power of two is exact. A norm that is zero, */ +/* infinite or NaN takes no part. DLASCL is called with both */ +/* endpoints at or above one so that forming the factor raises no */ +/* underflow. */ + + anrm = dlange_("M", m, n, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0. && anrm <= myhuge_(&c_b24)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b28, &i__3); + dlascl_("G", &c__0, &c__0, &r__1, &r__2, m, n, &a[a_offset], lda, + info); + } + + bnrm = dlange_("M", p, n, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0. && bnrm <= myhuge_(&c_b24)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b28, &i__3); + dlascl_("G", &c__0, &c__0, &r__1, &r__2, p, n, &b[b_offset], ldb, + info); + } + + i__1 = f2cmax(1,*m); + cnrm = dlange_("M", m, &c__1, &c__[1], &i__1, rwork); + i__1 = f2cmax(1,*p); + dnrm = dlange_("M", p, &c__1, &d__[1], &i__1, rwork); + ic = -myhuge_(&c__0); + if (cnrm > 0. && cnrm <= myhuge_(&c_b24)) { + ic = myexp_(&cnrm) + ka; + } + id = -myhuge_(&c__0); + if (dnrm > 0. && dnrm <= myhuge_(&c_b24)) { + id = myexp_(&dnrm) + kb; + } + iu = f2cmax(ic,id); + kt = 0; + if (iu > ibig) { + kt = ibig - iu; + } else if (iu < isml && (real) iu != -myhuge_(&c__0)) { + kt = isml - iu; + } + if (ka + kt != 0) { + i__1 = ka + kt; + r__1 = myscale_(&c_b28, &i__1); + dscal_(m, &r__1, &c__[1], &c__1); + } + if (kb + kt != 0) { + i__1 = kb + kt; + r__1 = myscale_(&c_b28, &i__1); + dscal_(p, &r__1, &d__[1], &c__1); + } + /* Compute the GRQ factorization of matrices B and A: */ /* B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P */ @@ -575,8 +683,8 @@ f"> */ /* Update c1 */ i__1 = *n - *p; - dgemv_("No transpose", &i__1, p, &c_b31, &a[(*n - *p + 1) * a_dim1 + - 1], lda, &d__[1], &c__1, &c_b33, &c__[1], &c__1); + dgemv_("No transpose", &i__1, p, &c_b60, &a[(*n - *p + 1) * a_dim1 + + 1], lda, &d__[1], &c__1, &c_b28, &c__[1], &c__1); } /* Solve R11*x1 = c1 for x1 */ @@ -604,8 +712,8 @@ f"> */ nr = *m + *p - *n; if (nr > 0) { i__1 = *n - *m; - dgemv_("No transpose", &nr, &i__1, &c_b31, &a[*n - *p + 1 + (*m + - 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b33, &c__[*n - + dgemv_("No transpose", &nr, &i__1, &c_b60, &a[*n - *p + 1 + (*m + + 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b28, &c__[*n - *p + 1], &c__1); } } else { @@ -614,7 +722,7 @@ f"> */ if (nr > 0) { dtrmv_("Upper", "No transpose", "Non unit", &nr, &a[*n - *p + 1 + (*n - *p + 1) * a_dim1], lda, &d__[1], &c__1); - daxpy_(&nr, &c_b31, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); + daxpy_(&nr, &c_b60, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); } /* Backward transformation x = Q**T*x */ @@ -622,6 +730,22 @@ f"> */ i__1 = *lwork - *p - mn; dormrq_("Left", "Transpose", n, &c__1, p, &b[b_offset], ldb, &work[1], &x[ 1], n, &work[*p + mn + 1], &i__1, info); + +/* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) */ +/* carries 2**(KA+KT) */ + + if (kt != 0) { + i__1 = -kt; + r__1 = myscale_(&c_b28, &i__1); + dscal_(n, &r__1, &x[1], &c__1); + } + if (ka + kt != 0 && *m + *p > *n) { + i__1 = *m + *p - *n; + i__2 = -ka - kt; + r__1 = myscale_(&c_b28, &i__2); + dscal_(&i__1, &r__1, &c__[*n - *p + 1], &c__1); + } + /* Computing MAX */ i__1 = lopt, i__2 = (integer) work[*p + mn + 1]; work[1] = (doublereal) (*p + mn + f2cmax(i__1,i__2)); diff --git a/lapack-netlib/SRC/dgglse.f b/lapack-netlib/SRC/dgglse.f index 28aeaf6e76..23db6b0c5a 100644 --- a/lapack-netlib/SRC/dgglse.f +++ b/lapack-netlib/SRC/dgglse.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download DGGLSE + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -52,6 +50,16 @@ *> matrices (B, A) given by *> *> B = (0 R)*Q, A = Z*T*Q. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The DTRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized RQ +*> factorization of the pair (B, A) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -153,12 +161,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with B in the -*> generalized RQ factorization of the pair (B, A) is +*> generalized RQ factorization of the pair (B, A) is exactly *> singular, so that rank(B) < P; the least squares *> solution could not be computed. *> = 2: the (N-P) by (N-P) part of the upper trapezoidal factor *> T associated with A in the generalized RQ factorization -*> of the pair (B, A) is singular, so that +*> of the pair (B, A) is exactly singular, so that *> rank( (A) ) < N; the least squares solution could not *> ( (B) ) *> be computed. @@ -172,11 +180,13 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup doubleOTHERsolve +*> \ingroup gglse * * ===================================================================== - SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, + SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -193,24 +203,30 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * ===================================================================== * * .. Parameters .. - DOUBLE PRECISION ONE - PARAMETER ( ONE = 1.0D+0 ) + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, $ NB4, NR + INTEGER IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT + DOUBLE PRECISION ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL DAXPY, DCOPY, DGEMV, DGGRQF, DORMQR, DORMRQ, - $ DTRMV, DTRTRS, XERBLA + EXTERNAL DAXPY, DCOPY, DGEMV, DGGRQF, DLASCL, DORMQR, + $ DORMRQ, DSCAL, DTRMV, DTRTRS, XERBLA * .. * .. External Functions .. INTEGER ILAENV - EXTERNAL ILAENV + DOUBLE PRECISION DLAMCH, DLANGE + EXTERNAL DLAMCH, DLANGE, ILAENV * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -265,6 +281,72 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( N.EQ.0 ) $ RETURN * +* Get machine parameters +* + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B, c and d by powers of two so that their largest +* entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies +* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and +* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which +* keeps the constraint and the residual consistent and scales x by +* 2**KT; KT brings the larger of the two right-hand sides into +* range. Scaling by a power of two is exact. A norm that is zero, +* infinite or NaN takes no part. DLASCL is called with both +* endpoints at or above one so that forming the factor raises no +* underflow. +* + ANRM = DLANGE( 'M', M, N, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO ) +* + BNRM = DLANGE( 'M', P, N, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL DLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO ) +* + CNRM = DLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK ) + DNRM = DLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK ) + IC = -HUGE( 0 ) + IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) ) + $ IC = EXPONENT( CNRM ) + KA + ID = -HUGE( 0 ) + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) + $ ID = EXPONENT( DNRM ) + KB + IU = MAX( IC, ID ) + KT = 0 + IF( IU.GT.IBIG ) THEN + KT = IBIG - IU + ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN + KT = ISML - IU + END IF + IF( KA+KT.NE.0 ) + $ CALL DSCAL( M, SCALE( ONE, KA+KT ), C, 1 ) + IF( KB+KT.NE.0 ) + $ CALL DSCAL( P, SCALE( ONE, KB+KT ), D, 1 ) +* * Compute the GRQ factorization of matrices B and A: * * B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P @@ -281,7 +363,8 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * Update c = Z**T *c = ( c1 ) N-P * ( c2 ) M+P-N * - CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ), + CALL DORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, + $ WORK( P+1 ), $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO ) LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) ) * @@ -302,7 +385,8 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Update c1 * - CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA, + CALL DGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), + $ LDA, $ D, 1, ONE, C, 1 ) END IF * @@ -327,7 +411,8 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( M.LT.N ) THEN NR = M + P - N IF( NR.GT.0 ) - $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ), + $ CALL DGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, + $ M+1 ), $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 ) ELSE NR = P @@ -340,8 +425,18 @@ SUBROUTINE DGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Backward transformation x = Q**T*x * - CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X, + CALL DORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), + $ X, $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) +* +* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) +* carries 2**(KA+KT) +* + IF( KT.NE.0 ) + $ CALL DSCAL( N, SCALE( ONE, -KT ), X, 1 ) + IF( KA+KT.NE.0 .AND. M+P.GT.N ) + $ CALL DSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 ) +* WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) * RETURN diff --git a/lapack-netlib/SRC/sggglm.c b/lapack-netlib/SRC/sggglm.c index c122bb2bff..2b5002c058 100644 --- a/lapack-netlib/SRC/sggglm.c +++ b/lapack-netlib/SRC/sggglm.c @@ -212,7 +212,10 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(float *x) {int e; (void)frexp(*x,&e); return e;} +static float myscale_(float *w, integer *i) {float myscale=*w*pow(2.,*(double*)i); return(myscale);} +static float myhuge_(float *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -222,14 +225,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c_n1 = -1; -static real c_b32 = -1.f; -static real c_b34 = 1.f; +static real c_b24 = 0.f; +static integer c__0 = 0; +static real c_b28 = 1.f; +static real c_b53 = -1.f; /* > \brief \b SGGGLM */ @@ -238,7 +241,6 @@ static real c_b34 = 1.f; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download SGGGLM + dependencies */ /* > */ @@ -249,7 +251,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -293,6 +294,16 @@ f"> */ /* > x */ /* > */ /* > where inv(B) denotes the inverse of B. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The STRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized QR */ +/* > factorization of the pair (A, B) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -391,12 +402,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with A in the */ -/* > generalized QR factorization of the pair (A, B) is */ +/* > generalized QR factorization of the pair (A, B) is exactly */ /* > singular, so that rank(A) < M; the least squares */ /* > solution could not be computed. */ /* > = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal */ /* > factor T associated with B in the generalized QR */ -/* > factorization of the pair (A, B) is singular, so that */ +/* > factorization of the pair (A, B) is exactly singular, so that */ /* > rank( A B ) < N; the least squares solution could not */ /* > be computed. */ /* > \endverbatim */ @@ -409,9 +420,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup realOTHEReigen */ +/* > \ingroup ggglm */ /* ===================================================================== */ /* Subroutine */ void sggglm_(integer *n, integer *m, integer *p, real *a, @@ -420,33 +429,44 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; + real r__1, r__2; /* Local variables */ - integer lopt, i__; - extern /* Subroutine */ void sgemv_(char *, integer *, integer *, real *, - real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *); - integer nb, np; + integer i__, ia, ib, ka, id, kb, kd, nb, np, nb1, nb2, nb3, nb4, ibig; + real anrm, bnrm, dnrm; + integer isml, lopt; + extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *), + sgemv_(char *, integer *, integer *, real *, real *, integer *, + real *, integer *, real *, real *, integer *), scopy_( + integer *, real *, integer *, real *, integer *); + real rwork[1]; + extern real slamch_(char *), slange_(char *, integer *, integer *, + real *, integer *, real *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - extern /* Subroutine */ void sggqrf_(integer *, integer *, integer *, real - *, integer *, real *, real *, integer *, real *, real *, integer * - , integer *); - integer lwkmin, nb1, nb2, nb3, nb4, lwkopt; + real bignum; + extern /* Subroutine */ void slascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, real *, integer *, integer *), sggqrf_(integer *, integer *, integer *, real *, integer + *, real *, real *, integer *, real *, real *, integer *, integer * + ); + integer lwkmin; + real smlnum; + integer lwkopt; logical lquery; extern /* Subroutine */ void sormqr_(char *, char *, integer *, integer *, integer *, real *, integer *, real *, real *, integer *, real *, integer *, integer *), sormrq_(char *, char *, integer *, integer *, integer *, real *, integer *, real *, real * - , integer *, real *, integer *, integer *); - extern void strtrs_(char *, char *, char *, integer *, integer *, real *, + , integer *, real *, integer *, integer *), + strtrs_(char *, char *, char *, integer *, integer *, real *, integer *, real *, integer *, integer *); + extern real sroundup_lwork__(integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* =================================================================== */ @@ -456,10 +476,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --d__; --x; @@ -503,7 +523,7 @@ f"> */ lwkmin = *m + *n + *p; lwkopt = *m + np + f2cmax(*n,*p) * nb; } - work[1] = (real) lwkopt; + work[1] = sroundup_lwork__(&lwkopt); if (*lwork < lwkmin && ! lquery) { *info = -12; @@ -532,6 +552,78 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = slamch_("S") / slamch_("P"); + bignum = 1.f / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B and d independently by powers of two 2**KA, 2**KB and */ +/* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the */ +/* values whose MYEXP lies in [ISML,IBIG]. x is then scaled by */ +/* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is */ +/* exact. A norm that is zero, infinite or NaN takes no part. */ +/* SLASCL is called with both endpoints at or above one so that */ +/* forming the factor raises no underflow. */ + + anrm = slange_("M", n, m, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0.f && anrm <= myhuge_(&c_b24)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b28, &i__3); + slascl_("G", &c__0, &c__0, &r__1, &r__2, n, m, &a[a_offset], lda, + info); + } + + bnrm = slange_("M", n, p, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0.f && bnrm <= myhuge_(&c_b24)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b28, &i__3); + slascl_("G", &c__0, &c__0, &r__1, &r__2, n, p, &b[b_offset], ldb, + info); + } + + dnrm = slange_("M", n, &c__1, &d__[1], n, rwork); + kd = 0; + if (dnrm > 0.f && dnrm <= myhuge_(&c_b24)) { + id = myexp_(&dnrm); + if (id < isml) { + kd = isml - id; + } else if (id > ibig) { + kd = ibig - id; + } + } + if (kd != 0) { + r__1 = myscale_(&c_b28, &kd); + sscal_(n, &r__1, &d__[1], &c__1); + } + /* Compute the GQR factorization of matrices A and B: */ /* Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M */ @@ -544,7 +636,7 @@ f"> */ i__1 = *lwork - *m - np; sggqrf_(n, m, p, &a[a_offset], lda, &work[1], &b[b_offset], ldb, &work[*m + 1], &work[*m + np + 1], &i__1, info); - lopt = work[*m + np + 1]; + lopt = (integer) work[*m + np + 1]; /* Update left-hand-side vector d = Q**T*d = ( d1 ) M */ /* ( d2 ) N-M */ @@ -586,8 +678,8 @@ f"> */ /* Update d1 = d1 - T12*y2 */ i__1 = *n - *m; - sgemv_("No transpose", m, &i__1, &c_b32, &b[(*m + *p - *n + 1) * b_dim1 + - 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b34, &d__[1], &c__1); + sgemv_("No transpose", m, &i__1, &c_b53, &b[(*m + *p - *n + 1) * b_dim1 + + 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b28, &d__[1], &c__1); /* Solve triangular system: R11*x = d1 */ @@ -613,6 +705,20 @@ f"> */ i__4 = *lwork - *m - np; sormrq_("Left", "Transpose", p, &c__1, &np, &b[f2cmax(i__1,i__2) + b_dim1], ldb, &work[*m + 1], &y[1], &i__3, &work[*m + np + 1], &i__4, info); + +/* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) */ + + if (ka != kd) { + i__1 = ka - kd; + r__1 = myscale_(&c_b28, &i__1); + sscal_(m, &r__1, &x[1], &c__1); + } + if (kb != kd) { + i__1 = kb - kd; + r__1 = myscale_(&c_b28, &i__1); + sscal_(p, &r__1, &y[1], &c__1); + } + /* Computing MAX */ i__1 = lopt, i__2 = (integer) work[*m + np + 1]; work[1] = (real) (*m + np + f2cmax(i__1,i__2)); diff --git a/lapack-netlib/SRC/sggglm.f b/lapack-netlib/SRC/sggglm.f index 37094e4f26..caed7cf5d1 100644 --- a/lapack-netlib/SRC/sggglm.f +++ b/lapack-netlib/SRC/sggglm.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SGGGLM + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -61,6 +59,16 @@ *> x *> *> where inv(B) denotes the inverse of B. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The STRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized QR +*> factorization of the pair (A, B) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -159,12 +167,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with A in the -*> generalized QR factorization of the pair (A, B) is +*> generalized QR factorization of the pair (A, B) is exactly *> singular, so that rank(A) < M; the least squares *> solution could not be computed. *> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal *> factor T associated with B in the generalized QR -*> factorization of the pair (A, B) is singular, so that +*> factorization of the pair (A, B) is exactly singular, so that *> rank( A B ) < N; the least squares solution could not *> be computed. *> \endverbatim @@ -180,8 +188,10 @@ *> \ingroup ggglm * * ===================================================================== - SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, + SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -205,18 +215,24 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, LOGICAL LQUERY INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, $ NB4, NP + INTEGER IA, IB, IBIG, ID, ISML, KA, KB, KD + REAL ANRM, BIGNUM, BNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + REAL RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL SCOPY, SGEMV, SGGQRF, SORMQR, SORMRQ, STRTRS, - $ XERBLA + EXTERNAL SCOPY, SGEMV, SGGQRF, SLASCL, SORMQR, SORMRQ, + $ SSCAL, STRTRS, XERBLA * .. * .. External Functions .. INTEGER ILAENV REAL SROUNDUP_LWORK - EXTERNAL ILAENV, SROUNDUP_LWORK + REAL SLAMCH, SLANGE + EXTERNAL ILAENV, SLAMCH, SLANGE, SROUNDUP_LWORK * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -278,6 +294,62 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, RETURN END IF * +* Get machine parameters +* + SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B and d independently by powers of two 2**KA, 2**KB and +* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the +* values whose EXPONENT lies in [ISML,IBIG]. x is then scaled by +* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is +* exact. A norm that is zero, infinite or NaN takes no part. +* SLASCL is called with both endpoints at or above one so that +* forming the factor raises no underflow. +* + ANRM = SLANGE( 'M', N, M, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO ) +* + BNRM = SLANGE( 'M', N, P, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO ) +* + DNRM = SLANGE( 'M', N, 1, D, N, RWORK ) + KD = 0 + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN + ID = EXPONENT( DNRM ) + IF( ID.LT.ISML ) THEN + KD = ISML - ID + ELSE IF( ID.GT.IBIG ) THEN + KD = IBIG - ID + END IF + END IF + IF( KD.NE.0 ) + $ CALL SSCAL( N, SCALE( ONE, KD ), D, 1 ) +* * Compute the GQR factorization of matrices A and B: * * Q**T*A = ( R11 ) M, Q**T*B*Z**T = ( T11 T12 ) M @@ -326,7 +398,8 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, * Solve triangular system: R11*x = d1 * IF( M.GT.0 ) THEN - CALL STRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, + CALL STRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, + $ LDA, $ D, M, INFO ) * IF( INFO.GT.0 ) THEN @@ -344,7 +417,15 @@ SUBROUTINE SGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, CALL SORMRQ( 'Left', 'Transpose', P, 1, NP, $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) - WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) +* +* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) +* + IF( KA.NE.KD ) + $ CALL SSCAL( M, SCALE( ONE, KA-KD ), X, 1 ) + IF( KB.NE.KD ) + $ CALL SSCAL( P, SCALE( ONE, KB-KD ), Y, 1 ) +* + WORK( 1 ) = REAL( M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) ) * RETURN * diff --git a/lapack-netlib/SRC/sgglse.c b/lapack-netlib/SRC/sgglse.c index 75e379a427..74b0e0a161 100644 --- a/lapack-netlib/SRC/sgglse.c +++ b/lapack-netlib/SRC/sgglse.c @@ -212,7 +212,11 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(float *x) {int e; (void)frexp(*x,&e); return e;} +static float myscale_(float *w, integer *i) {float myscale=*w*pow(2.,*(double*)i); return(myscale);} +static float myhuge_(float *w) {return(HUGE_VAL);} +static int myihuge_(int *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -222,14 +226,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static integer c__1 = 1; static integer c_n1 = -1; -static real c_b31 = -1.f; -static real c_b33 = 1.f; +static real c_b24 = 0.f; +static integer c__0 = 0; +static real c_b28 = 1.f; +static real c_b60 = -1.f; /* > \brief SGGLSE solves overdetermined or underdetermined systems for OTHER matrices */ @@ -238,7 +242,6 @@ static real c_b33 = 1.f; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download SGGLSE + dependencies */ /* > */ @@ -249,7 +252,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -284,6 +286,16 @@ f"> */ /* > matrices (B, A) given by */ /* > */ /* > B = (0 R)*Q, A = Z*T*Q. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The STRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized RQ */ +/* > factorization of the pair (B, A) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -385,12 +397,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with B in the */ -/* > generalized RQ factorization of the pair (B, A) is */ +/* > generalized RQ factorization of the pair (B, A) is exactly */ /* > singular, so that rank(B) < P; the least squares */ /* > solution could not be computed. */ /* > = 2: the (N-P) by (N-P) part of the upper trapezoidal factor */ /* > T associated with A in the generalized RQ factorization */ -/* > of the pair (B, A) is singular, so that */ +/* > of the pair (B, A) is exactly singular, so that */ /* > rank( (A) ) < N; the least squares solution could not */ /* > ( (B) ) */ /* > be computed. */ @@ -404,9 +416,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup realOTHERsolve */ +/* > \ingroup gglse */ /* ===================================================================== */ /* Subroutine */ void sgglse_(integer *m, integer *n, integer *p, real *a, @@ -414,37 +424,49 @@ f"> */ real *work, integer *lwork, integer *info) { /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; + integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3; + real r__1, r__2; /* Local variables */ - integer lopt; - extern /* Subroutine */ void sgemv_(char *, integer *, integer *, real *, - real *, integer *, real *, integer *, real *, real *, integer *), scopy_(integer *, real *, integer *, real *, integer *), - saxpy_(integer *, real *, real *, integer *, real *, integer *), - strmv_(char *, char *, char *, integer *, real *, integer *, real - *, integer *); - integer nb, mn, nr; + integer ia, ib, ic, id, ka, kb, nb, mn, iu, kt, nr, nb1, nb2, nb3, nb4, + ibig; + real anrm, bnrm, cnrm, dnrm; + integer isml, lopt; + extern /* Subroutine */ void sscal_(integer *, real *, real *, integer *), + sgemv_(char *, integer *, integer *, real *, real *, integer *, + real *, integer *, real *, real *, integer *), scopy_( + integer *, real *, integer *, real *, integer *); + real rwork[1]; + extern /* Subroutine */ void saxpy_(integer *, real *, real *, integer *, + real *, integer *), strmv_(char *, char *, char *, integer *, + real *, integer *, real *, integer *); + extern real slamch_(char *), slange_(char *, integer *, integer *, + real *, integer *, real *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); - extern /* Subroutine */ void sggrqf_(integer *, integer *, integer *, real - *, integer *, real *, real *, integer *, real *, real *, integer * - , integer *); - integer lwkmin, nb1, nb2, nb3, nb4, lwkopt; + real bignum; + extern /* Subroutine */ void slascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, real *, integer *, integer *), sggrqf_(integer *, integer *, integer *, real *, integer + *, real *, real *, integer *, real *, real *, integer *, integer * + ); + integer lwkmin; + real smlnum; + integer lwkopt; logical lquery; extern /* Subroutine */ void sormqr_(char *, char *, integer *, integer *, integer *, real *, integer *, real *, real *, integer *, real *, integer *, integer *), sormrq_(char *, char *, integer *, integer *, integer *, real *, integer *, real *, real * - , integer *, real *, integer *, integer *); - extern void strtrs_(char *, char *, char *, integer *, integer *, real *, + , integer *, real *, integer *, integer *), + strtrs_(char *, char *, char *, integer *, integer *, real *, integer *, real *, integer *, integer *); + extern real sroundup_lwork__(integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* ===================================================================== */ @@ -454,10 +476,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --c__; --d__; @@ -501,7 +523,7 @@ f"> */ lwkmin = *m + *n + *p; lwkopt = *p + mn + f2cmax(*m,*n) * nb; } - work[1] = (real) lwkopt; + work[1] = sroundup_lwork__(&lwkopt); if (*lwork < lwkmin && ! lquery) { *info = -12; @@ -522,6 +544,96 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = slamch_("S") / slamch_("P"); + bignum = 1.f / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B, c and d by powers of two so that their largest */ +/* entries lie in [SMLNUM,BIGNUM), the values whose MYEXP lies */ +/* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and */ +/* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which */ +/* keeps the constraint and the residual consistent and scales x by */ +/* 2**KT; KT brings the larger of the two right-hand sides into */ +/* range. Scaling by a power of two is exact. A norm that is zero, */ +/* infinite or NaN takes no part. SLASCL is called with both */ +/* endpoints at or above one so that forming the factor raises no */ +/* underflow. */ + + anrm = slange_("M", m, n, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0.f && anrm <= myhuge_(&c_b24)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b28, &i__3); + slascl_("G", &c__0, &c__0, &r__1, &r__2, m, n, &a[a_offset], lda, + info); + } + + bnrm = slange_("M", p, n, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0.f && bnrm <= myhuge_(&c_b24)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b28, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b28, &i__3); + slascl_("G", &c__0, &c__0, &r__1, &r__2, p, n, &b[b_offset], ldb, + info); + } + + i__1 = f2cmax(1,*m); + cnrm = slange_("M", m, &c__1, &c__[1], &i__1, rwork); + i__1 = f2cmax(1,*p); + dnrm = slange_("M", p, &c__1, &d__[1], &i__1, rwork); + ic = -myihuge_(&c__0); + if (cnrm > 0.f && cnrm <= myhuge_(&c_b24)) { + ic = myexp_(&cnrm) + ka; + } + id = -myihuge_(&c__0); + if (dnrm > 0.f && dnrm <= myhuge_(&c_b24)) { + id = myexp_(&dnrm) + kb; + } + iu = f2cmax(ic,id); + kt = 0; + if (iu > ibig) { + kt = ibig - iu; + } else if (iu < isml && (real) iu != -myihuge_(&c__0)) { + kt = isml - iu; + } + if (ka + kt != 0) { + i__1 = ka + kt; + r__1 = myscale_(&c_b28, &i__1); + sscal_(m, &r__1, &c__[1], &c__1); + } + if (kb + kt != 0) { + i__1 = kb + kt; + r__1 = myscale_(&c_b28, &i__1); + sscal_(p, &r__1, &d__[1], &c__1); + } + /* Compute the GRQ factorization of matrices B and A: */ /* B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P */ @@ -534,7 +646,7 @@ f"> */ i__1 = *lwork - *p - mn; sggrqf_(p, m, n, &b[b_offset], ldb, &work[1], &a[a_offset], lda, &work[*p + 1], &work[*p + mn + 1], &i__1, info); - lopt = work[*p + mn + 1]; + lopt = (integer) work[*p + mn + 1]; /* Update c = Z**T *c = ( c1 ) N-P */ /* ( c2 ) M+P-N */ @@ -565,8 +677,8 @@ f"> */ /* Update c1 */ i__1 = *n - *p; - sgemv_("No transpose", &i__1, p, &c_b31, &a[(*n - *p + 1) * a_dim1 + - 1], lda, &d__[1], &c__1, &c_b33, &c__[1], &c__1); + sgemv_("No transpose", &i__1, p, &c_b60, &a[(*n - *p + 1) * a_dim1 + + 1], lda, &d__[1], &c__1, &c_b28, &c__[1], &c__1); } /* Solve R11*x1 = c1 for x1 */ @@ -594,8 +706,8 @@ f"> */ nr = *m + *p - *n; if (nr > 0) { i__1 = *n - *m; - sgemv_("No transpose", &nr, &i__1, &c_b31, &a[*n - *p + 1 + (*m + - 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b33, &c__[*n - + sgemv_("No transpose", &nr, &i__1, &c_b60, &a[*n - *p + 1 + (*m + + 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b28, &c__[*n - *p + 1], &c__1); } } else { @@ -604,7 +716,7 @@ f"> */ if (nr > 0) { strmv_("Upper", "No transpose", "Non unit", &nr, &a[*n - *p + 1 + (*n - *p + 1) * a_dim1], lda, &d__[1], &c__1); - saxpy_(&nr, &c_b31, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); + saxpy_(&nr, &c_b60, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); } /* Backward transformation x = Q**T*x */ @@ -612,6 +724,22 @@ f"> */ i__1 = *lwork - *p - mn; sormrq_("Left", "Transpose", n, &c__1, p, &b[b_offset], ldb, &work[1], &x[ 1], n, &work[*p + mn + 1], &i__1, info); + +/* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) */ +/* carries 2**(KA+KT) */ + + if (kt != 0) { + i__1 = -kt; + r__1 = myscale_(&c_b28, &i__1); + sscal_(n, &r__1, &x[1], &c__1); + } + if (ka + kt != 0 && *m + *p > *n) { + i__1 = *m + *p - *n; + i__2 = -ka - kt; + r__1 = myscale_(&c_b28, &i__2); + sscal_(&i__1, &r__1, &c__[*n - *p + 1], &c__1); + } + /* Computing MAX */ i__1 = lopt, i__2 = (integer) work[*p + mn + 1]; work[1] = (real) (*p + mn + f2cmax(i__1,i__2)); diff --git a/lapack-netlib/SRC/sgglse.f b/lapack-netlib/SRC/sgglse.f index 53e3f8e45b..1f2338c579 100644 --- a/lapack-netlib/SRC/sgglse.f +++ b/lapack-netlib/SRC/sgglse.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download SGGLSE + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -52,6 +50,16 @@ *> matrices (B, A) given by *> *> B = (0 R)*Q, A = Z*T*Q. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The STRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized RQ +*> factorization of the pair (B, A) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -153,12 +161,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with B in the -*> generalized RQ factorization of the pair (B, A) is +*> generalized RQ factorization of the pair (B, A) is exactly *> singular, so that rank(B) < P; the least squares *> solution could not be computed. *> = 2: the (N-P) by (N-P) part of the upper trapezoidal factor *> T associated with A in the generalized RQ factorization -*> of the pair (B, A) is singular, so that +*> of the pair (B, A) is exactly singular, so that *> rank( (A) ) < N; the least squares solution could not *> ( (B) ) *> be computed. @@ -175,8 +183,10 @@ *> \ingroup gglse * * ===================================================================== - SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, + SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -193,25 +203,31 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * ===================================================================== * * .. Parameters .. - REAL ONE - PARAMETER ( ONE = 1.0E+0 ) + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) * .. * .. Local Scalars .. LOGICAL LQUERY INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, $ NB4, NR + INTEGER IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT + REAL ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + REAL RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL SAXPY, SCOPY, SGEMV, SGGRQF, SORMQR, SORMRQ, - $ STRMV, STRTRS, XERBLA + EXTERNAL SAXPY, SCOPY, SGEMV, SGGRQF, SLASCL, SORMQR, + $ SORMRQ, SSCAL, STRMV, STRTRS, XERBLA * .. * .. External Functions .. INTEGER ILAENV REAL SROUNDUP_LWORK - EXTERNAL ILAENV, SROUNDUP_LWORK + REAL SLAMCH, SLANGE + EXTERNAL ILAENV, SLAMCH, SLANGE, SROUNDUP_LWORK * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -266,6 +282,72 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( N.EQ.0 ) $ RETURN * +* Get machine parameters +* + SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B, c and d by powers of two so that their largest +* entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies +* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and +* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which +* keeps the constraint and the residual consistent and scales x by +* 2**KT; KT brings the larger of the two right-hand sides into +* range. Scaling by a power of two is exact. A norm that is zero, +* infinite or NaN takes no part. SLASCL is called with both +* endpoints at or above one so that forming the factor raises no +* underflow. +* + ANRM = SLANGE( 'M', M, N, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO ) +* + BNRM = SLANGE( 'M', P, N, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL SLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO ) +* + CNRM = SLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK ) + DNRM = SLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK ) + IC = -HUGE( 0 ) + IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) ) + $ IC = EXPONENT( CNRM ) + KA + ID = -HUGE( 0 ) + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) + $ ID = EXPONENT( DNRM ) + KB + IU = MAX( IC, ID ) + KT = 0 + IF( IU.GT.IBIG ) THEN + KT = IBIG - IU + ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN + KT = ISML - IU + END IF + IF( KA+KT.NE.0 ) + $ CALL SSCAL( M, SCALE( ONE, KA+KT ), C, 1 ) + IF( KB+KT.NE.0 ) + $ CALL SSCAL( P, SCALE( ONE, KB+KT ), D, 1 ) +* * Compute the GRQ factorization of matrices B and A: * * B*Q**T = ( 0 T12 ) P Z**T*A*Q**T = ( R11 R12 ) N-P @@ -282,7 +364,8 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * Update c = Z**T *c = ( c1 ) N-P * ( c2 ) M+P-N * - CALL SORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, WORK( P+1 ), + CALL SORMQR( 'Left', 'Transpose', M, 1, MN, A, LDA, + $ WORK( P+1 ), $ C, MAX( 1, M ), WORK( P+MN+1 ), LWORK-P-MN, INFO ) LOPT = MAX( LOPT, INT( WORK( P+MN+1 ) ) ) * @@ -303,7 +386,8 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Update c1 * - CALL SGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), LDA, + CALL SGEMV( 'No transpose', N-P, P, -ONE, A( 1, N-P+1 ), + $ LDA, $ D, 1, ONE, C, 1 ) END IF * @@ -328,7 +412,8 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( M.LT.N ) THEN NR = M + P - N IF( NR.GT.0 ) - $ CALL SGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, M+1 ), + $ CALL SGEMV( 'No transpose', NR, N-M, -ONE, A( N-P+1, + $ M+1 ), $ LDA, D( NR+1 ), 1, ONE, C( N-P+1 ), 1 ) ELSE NR = P @@ -341,9 +426,19 @@ SUBROUTINE SGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Backward transformation x = Q**T*x * - CALL SORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), X, + CALL SORMRQ( 'Left', 'Transpose', N, 1, P, B, LDB, WORK( 1 ), + $ X, $ N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) - WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) +* +* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) +* carries 2**(KA+KT) +* + IF( KT.NE.0 ) + $ CALL SSCAL( N, SCALE( ONE, -KT ), X, 1 ) + IF( KA+KT.NE.0 .AND. M+P.GT.N ) + $ CALL SSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 ) +* + WORK( 1 ) = REAL( P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) ) * RETURN * diff --git a/lapack-netlib/SRC/zggglm.c b/lapack-netlib/SRC/zggglm.c index 06aeb31c12..dfde9dd23a 100644 --- a/lapack-netlib/SRC/zggglm.c +++ b/lapack-netlib/SRC/zggglm.c @@ -240,7 +240,10 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(double *x) {int e; (void)frexp(*x,&e); return e;} +static double myscale_(double *w, integer *i) {double myscale=*w*pow(2.,*(double*)i); return(myscale);} +static double myhuge_(double *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -250,13 +253,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static doublecomplex c_b2 = {1.,0.}; static integer c__1 = 1; static integer c_n1 = -1; +static doublereal c_b26 = 0.; +static integer c__0 = 0; +static doublereal c_b30 = 1.; /* > \brief \b ZGGGLM */ @@ -265,7 +269,6 @@ static integer c_n1 = -1; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download ZGGGLM + dependencies */ /* > */ @@ -276,7 +279,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -320,6 +322,16 @@ f"> */ /* > x */ /* > */ /* > where inv(B) denotes the inverse of B. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The ZTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized QR */ +/* > factorization of the pair (A, B) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -418,12 +430,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with A in the */ -/* > generalized QR factorization of the pair (A, B) is */ +/* > generalized QR factorization of the pair (A, B) is exactly */ /* > singular, so that rank(A) < M; the least squares */ /* > solution could not be computed. */ /* > = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal */ /* > factor T associated with B in the generalized QR */ -/* > factorization of the pair (A, B) is singular, so that */ +/* > factorization of the pair (A, B) is exactly singular, so that */ /* > rank( A B ) < N; the least squares solution could not */ /* > be computed. */ /* > \endverbatim */ @@ -436,9 +448,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complex16OTHEReigen */ +/* > \ingroup ggglm */ /* ===================================================================== */ /* Subroutine */ void zggglm_(integer *n, integer *m, integer *p, @@ -448,39 +458,51 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; + real r__1, r__2; doublecomplex z__1; /* Local variables */ - integer lopt, i__; + integer i__, ia, ib, ka, id, kb, kd, nb, np, nb1, nb2, nb3, nb4, ibig; + doublereal anrm, bnrm, dnrm; + integer isml, lopt; extern /* Subroutine */ void zgemv_(char *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, - integer *, doublecomplex *, doublecomplex *, integer *), - zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, - integer *); - integer nb, np; + integer *, doublecomplex *, doublecomplex *, integer *); + doublereal rwork[1]; + extern /* Subroutine */ void zcopy_(integer *, doublecomplex *, integer *, + doublecomplex *, integer *); + extern doublereal dlamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); + extern /* Subroutine */ void zdscal_(integer *, real *, doublecomplex *, + integer *); + doublereal bignum; + extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, + integer *, doublereal *); + extern /* Subroutine */ void zlascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, doublecomplex *, integer *, integer + *); extern /* Subroutine */ void zggqrf_(integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, integer *) ; - integer lwkmin, nb1, nb2, nb3, nb4, lwkopt; + integer lwkmin; + doublereal smlnum; + integer lwkopt; logical lquery; extern /* Subroutine */ void zunmqr_(char *, char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, integer *, integer *), zunmrq_(char *, char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, - doublecomplex *, integer *, doublecomplex *, integer *, integer *); - extern void ztrtrs_(char *, char *, char *, integer *, + doublecomplex *, integer *, doublecomplex *, integer *, integer *), ztrtrs_(char *, char *, char *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, integer *, integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* =================================================================== */ @@ -490,10 +512,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --d__; --x; @@ -568,6 +590,78 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = dlamch_("S") / dlamch_("P"); + bignum = 1. / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B and d independently by powers of two 2**KA, 2**KB and */ +/* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the */ +/* values whose MYEXP lies in [ISML,IBIG]. x is then scaled by */ +/* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is */ +/* exact. A norm that is zero, infinite or NaN takes no part. */ +/* ZLASCL is called with both endpoints at or above one so that */ +/* forming the factor raises no underflow. */ + + anrm = zlange_("M", n, m, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0. && anrm <= myhuge_(&c_b26)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b30, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b30, &i__3); + zlascl_("G", &c__0, &c__0, &r__1, &r__2, n, m, &a[a_offset], lda, + info); + } + + bnrm = zlange_("M", n, p, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0. && bnrm <= myhuge_(&c_b26)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b30, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b30, &i__3); + zlascl_("G", &c__0, &c__0, &r__1, &r__2, n, p, &b[b_offset], ldb, + info); + } + + dnrm = zlange_("M", n, &c__1, &d__[1], n, rwork); + kd = 0; + if (dnrm > 0. && dnrm <= myhuge_(&c_b26)) { + id = myexp_(&dnrm); + if (id < isml) { + kd = isml - id; + } else if (id > ibig) { + kd = ibig - id; + } + } + if (kd != 0) { + r__1 = myscale_(&c_b30, &kd); + zdscal_(n, &r__1, &d__[1], &c__1); + } + /* Compute the GQR factorization of matrices A and B: */ /* Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M */ @@ -625,7 +719,7 @@ f"> */ /* Update d1 = d1 - T12*y2 */ i__1 = *n - *m; - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; zgemv_("No transpose", m, &i__1, &z__1, &b[(*m + *p - *n + 1) * b_dim1 + 1], ldb, &y[*m + *p - *n + 1], &c__1, &c_b2, &d__[1], &c__1); @@ -654,6 +748,20 @@ f"> */ zunmrq_("Left", "Conjugate transpose", p, &c__1, &np, &b[f2cmax(i__1,i__2) + b_dim1], ldb, &work[*m + 1], &y[1], &i__3, &work[*m + np + 1], & i__4, info); + +/* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) */ + + if (ka != kd) { + i__1 = ka - kd; + r__1 = myscale_(&c_b30, &i__1); + zdscal_(m, &r__1, &x[1], &c__1); + } + if (kb != kd) { + i__1 = kb - kd; + r__1 = myscale_(&c_b30, &i__1); + zdscal_(p, &r__1, &y[1], &c__1); + } + /* Computing MAX */ i__4 = *m + np + 1; i__2 = lopt, i__3 = (integer) work[i__4].r; diff --git a/lapack-netlib/SRC/zggglm.f b/lapack-netlib/SRC/zggglm.f index 62b4acdec3..6bced07488 100644 --- a/lapack-netlib/SRC/zggglm.f +++ b/lapack-netlib/SRC/zggglm.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZGGGLM + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -61,6 +59,16 @@ *> x *> *> where inv(B) denotes the inverse of B. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The ZTRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized QR +*> factorization of the pair (A, B) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -159,12 +167,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with A in the -*> generalized QR factorization of the pair (A, B) is +*> generalized QR factorization of the pair (A, B) is exactly *> singular, so that rank(A) < M; the least squares *> solution could not be computed. *> = 2: the bottom (N-M) by (N-M) part of the upper trapezoidal *> factor T associated with B in the generalized QR -*> factorization of the pair (A, B) is singular, so that +*> factorization of the pair (A, B) is exactly singular, so that *> rank( A B ) < N; the least squares solution could not *> be computed. *> \endverbatim @@ -177,11 +185,13 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHEReigen +*> \ingroup ggglm * * ===================================================================== - SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, + SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -198,6 +208,8 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, * =================================================================== * * .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) COMPLEX*16 CZERO, CONE PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), $ CONE = ( 1.0D+0, 0.0D+0 ) ) @@ -206,17 +218,23 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, LOGICAL LQUERY INTEGER I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3, $ NB4, NP + INTEGER IA, IB, IBIG, ID, ISML, KA, KB, KD + DOUBLE PRECISION ANRM, BIGNUM, BNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZCOPY, ZGEMV, ZGGQRF, ZTRTRS, ZUNMQR, - $ ZUNMRQ + EXTERNAL XERBLA, ZCOPY, ZDSCAL, ZGEMV, ZGGQRF, ZLASCL, + $ ZTRTRS, ZUNMQR, ZUNMRQ * .. * .. External Functions .. INTEGER ILAENV - EXTERNAL ILAENV + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL DLAMCH, ILAENV, ZLANGE * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -278,6 +296,62 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, RETURN END IF * +* Get machine parameters +* + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B and d independently by powers of two 2**KA, 2**KB and +* 2**KD so that their largest entries lie in [SMLNUM,BIGNUM), the +* values whose EXPONENT lies in [ISML,IBIG]. x is then scaled by +* 2**(KD-KA) and y by 2**(KD-KB). Scaling by a power of two is +* exact. A norm that is zero, infinite or NaN takes no part. +* ZLASCL is called with both endpoints at or above one so that +* forming the factor raises no underflow. +* + ANRM = ZLANGE( 'M', N, M, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), N, M, A, LDA, INFO ) +* + BNRM = ZLANGE( 'M', N, P, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), N, P, B, LDB, INFO ) +* + DNRM = ZLANGE( 'M', N, 1, D, N, RWORK ) + KD = 0 + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) THEN + ID = EXPONENT( DNRM ) + IF( ID.LT.ISML ) THEN + KD = ISML - ID + ELSE IF( ID.GT.IBIG ) THEN + KD = IBIG - ID + END IF + END IF + IF( KD.NE.0 ) + $ CALL ZDSCAL( N, SCALE( ONE, KD ), D, 1 ) +* * Compute the GQR factorization of matrices A and B: * * Q**H*A = ( R11 ) M, Q**H*B*Z**H = ( T11 T12 ) M @@ -294,7 +368,8 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, * Update left-hand-side vector d = Q**H*d = ( d1 ) M * ( d2 ) N-M * - CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, WORK, + CALL ZUNMQR( 'Left', 'Conjugate transpose', N, 1, M, A, LDA, + $ WORK, $ D, MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) ) * @@ -320,13 +395,15 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, * * Update d1 = d1 - T12*y2 * - CALL ZGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), LDB, + CALL ZGEMV( 'No transpose', M, N-M, -CONE, B( 1, M+P-N+1 ), + $ LDB, $ Y( M+P-N+1 ), 1, CONE, D, 1 ) * * Solve triangular system: R11*x = d1 * IF( M.GT.0 ) THEN - CALL ZTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, LDA, + CALL ZTRTRS( 'Upper', 'No Transpose', 'Non unit', M, 1, A, + $ LDA, $ D, M, INFO ) * IF( INFO.GT.0 ) THEN @@ -344,6 +421,14 @@ SUBROUTINE ZGGGLM( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK, CALL ZUNMRQ( 'Left', 'Conjugate transpose', P, 1, NP, $ B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y, $ MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO ) +* +* Undo scaling: x carries 2**(KD-KA) and y carries 2**(KD-KB) +* + IF( KA.NE.KD ) + $ CALL ZDSCAL( M, SCALE( ONE, KA-KD ), X, 1 ) + IF( KB.NE.KD ) + $ CALL ZDSCAL( P, SCALE( ONE, KB-KD ), Y, 1 ) +* WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) ) * RETURN diff --git a/lapack-netlib/SRC/zgglse.c b/lapack-netlib/SRC/zgglse.c index fc709cb55f..74bc3dfb0d 100644 --- a/lapack-netlib/SRC/zgglse.c +++ b/lapack-netlib/SRC/zgglse.c @@ -240,7 +240,10 @@ typedef struct Namelist Namelist; #define myexit_() break; #define mycycle() continue; #define myceiling(w) {ceil(w)} -#define myhuge(w) {HUGE_VAL} +static int my_expfunc(double *x) {int e; (void)frexp(*x,&e); return e;} +static double myscale_(double *w, integer *i) {double myscale=*w*pow(2.,*(double*)i); return(myscale);} +static double myhuge_(double *w) {return(HUGE_VAL);} + //#define mymaxloc_(w,s,e,n) {if (sizeof(*(w)) == sizeof(double)) dmaxloc_((w),*(s),*(e),n); else dmaxloc_((w),*(s),*(e),n);} #define mymaxloc(w,s,e,n) {dmaxloc_(w,*(s),*(e),n)} @@ -250,13 +253,14 @@ typedef struct Namelist Namelist; */ - - /* Table of constant values */ static doublecomplex c_b1 = {1.,0.}; static integer c__1 = 1; static integer c_n1 = -1; +static doublereal c_b25 = 0.; +static integer c__0 = 0; +static doublereal c_b29 = 1.; /* > \brief ZGGLSE solves overdetermined or underdetermined systems for OTHER matrices */ @@ -265,7 +269,6 @@ static integer c_n1 = -1; /* Online html documentation available at */ /* http://www.netlib.org/lapack/explore-html/ */ -/* > \htmlonly */ /* > Download ZGGLSE + dependencies */ /* > */ @@ -276,7 +279,6 @@ f"> */ /* > */ /* > [TXT] */ -/* > \endhtmlonly */ /* Definition: */ /* =========== */ @@ -311,6 +313,16 @@ f"> */ /* > matrices (B, A) given by */ /* > */ /* > B = (0 R)*Q, A = Z*T*Q. */ +/* > */ +/* > Callers of this subroutine should note that the singularity/rank-deficiency checks */ +/* > implemented in this subroutine are rudimentary. The ZTRTRS subroutine called by this */ +/* > subroutine only signals a failure due to singularity if the problem is exactly singular. */ +/* > */ +/* > It is conceivable for one (or more) of the factors involved in the generalized RQ */ +/* > factorization of the pair (B, A) to be subnormally close to singularity without this */ +/* > subroutine signalling an error. The solutions computed for such almost-rank-deficient */ +/* > problems may be less accurate due to a loss of numerical precision. */ +/* > */ /* > \endverbatim */ /* Arguments: */ @@ -412,12 +424,12 @@ f"> */ /* > = 0: successful exit. */ /* > < 0: if INFO = -i, the i-th argument had an illegal value. */ /* > = 1: the upper triangular factor R associated with B in the */ -/* > generalized RQ factorization of the pair (B, A) is */ +/* > generalized RQ factorization of the pair (B, A) is exactly */ /* > singular, so that rank(B) < P; the least squares */ /* > solution could not be computed. */ /* > = 2: the (N-P) by (N-P) part of the upper trapezoidal factor */ /* > T associated with A in the generalized RQ factorization */ -/* > of the pair (B, A) is singular, so that */ +/* > of the pair (B, A) is exactly singular, so that */ /* > rank( (A) ) < N; the least squares solution could not */ /* > ( (B) ) */ /* > be computed. */ @@ -431,9 +443,7 @@ f"> */ /* > \author Univ. of Colorado Denver */ /* > \author NAG Ltd. */ -/* > \date December 2016 */ - -/* > \ingroup complex16OTHERsolve */ +/* > \ingroup gglse */ /* ===================================================================== */ /* Subroutine */ void zgglse_(integer *m, integer *n, integer *p, @@ -443,27 +453,43 @@ f"> */ { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; + real r__1, r__2; doublecomplex z__1; /* Local variables */ - integer lopt; + integer ia, ib, ic, id, ka, kb, nb, mn, iu, kt, nr, nb1, nb2, nb3, nb4, + ibig; + doublereal anrm, bnrm, cnrm, dnrm; + integer isml, lopt; extern /* Subroutine */ void zgemv_(char *, integer *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, - integer *, doublecomplex *, doublecomplex *, integer *), - zcopy_(integer *, doublecomplex *, integer *, doublecomplex *, - integer *), zaxpy_(integer *, doublecomplex *, doublecomplex *, - integer *, doublecomplex *, integer *), ztrmv_(char *, char *, - char *, integer *, doublecomplex *, integer *, doublecomplex *, - integer *); - integer nb, mn, nr; + integer *, doublecomplex *, doublecomplex *, integer *); + extern integer myexp_(doublereal *); + doublereal rwork[1]; + extern /* Subroutine */ void zcopy_(integer *, doublecomplex *, integer *, + doublecomplex *, integer *), zaxpy_(integer *, doublecomplex *, + doublecomplex *, integer *, doublecomplex *, integer *), ztrmv_( + char *, char *, char *, integer *, doublecomplex *, integer *, + doublecomplex *, integer *); + extern doublereal dlamch_(char *); extern /* Subroutine */ void xerbla_(char *, integer *, ftnlen); extern integer ilaenv_(integer *, char *, char *, integer *, integer *, integer *, integer *, ftnlen, ftnlen); + extern /* Subroutine */ void zdscal_(integer *, real *, doublecomplex *, + integer *); + doublereal bignum; + extern doublereal zlange_(char *, integer *, integer *, doublecomplex *, + integer *, doublereal *); + extern /* Subroutine */ void zlascl_(char *, integer *, integer *, real *, + real *, integer *, integer *, doublecomplex *, integer *, integer + *); extern /* Subroutine */ void zggrqf_(integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, doublecomplex *, doublecomplex *, integer *, integer *) ; - integer lwkmin, nb1, nb2, nb3, nb4, lwkopt; + integer lwkmin; + doublereal smlnum; + integer lwkopt; logical lquery; extern /* Subroutine */ void zunmqr_(char *, char *, integer *, integer *, integer *, doublecomplex *, integer *, doublecomplex *, @@ -474,10 +500,9 @@ f"> */ integer *); -/* -- LAPACK driver routine (version 3.7.0) -- */ +/* -- LAPACK driver routine -- */ /* -- LAPACK is a software package provided by Univ. of Tennessee, -- */ /* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */ -/* December 2016 */ /* ===================================================================== */ @@ -487,10 +512,10 @@ f"> */ /* Parameter adjustments */ a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; + a_offset = 1 + a_dim1; a -= a_offset; b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; + b_offset = 1 + b_dim1; b -= b_offset; --c__; --d__; @@ -555,6 +580,96 @@ f"> */ return; } +/* Get machine parameters */ + + smlnum = dlamch_("S") / dlamch_("P"); + bignum = 1. / smlnum; + isml = myexp_(&smlnum); + ibig = myexp_(&bignum) - 1; + +/* Scale A, B, c and d by powers of two so that their largest */ +/* entries lie in [SMLNUM,BIGNUM), the values whose MYEXP lies */ +/* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and */ +/* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which */ +/* keeps the constraint and the residual consistent and scales x by */ +/* 2**KT; KT brings the larger of the two right-hand sides into */ +/* range. Scaling by a power of two is exact. A norm that is zero, */ +/* infinite or NaN takes no part. ZLASCL is called with both */ +/* endpoints at or above one so that forming the factor raises no */ +/* underflow. */ + + anrm = zlange_("M", m, n, &a[a_offset], lda, rwork); + ka = 0; + if (anrm > 0. && anrm <= myhuge_(&c_b25)) { + ia = myexp_(&anrm); + if (ia < isml) { + ka = isml - ia; + } else if (ia > ibig) { + ka = ibig - ia; + } + } + if (ka != 0) { +/* Computing MAX */ + i__2 = -ka; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b29, &i__1); + i__3 = f2cmax(ka,0); + r__2 = myscale_(&c_b29, &i__3); + zlascl_("G", &c__0, &c__0, &r__1, &r__2, m, n, &a[a_offset], lda, + info); + } + + bnrm = zlange_("M", p, n, &b[b_offset], ldb, rwork); + kb = 0; + if (bnrm > 0. && bnrm <= myhuge_(&c_b25)) { + ib = myexp_(&bnrm); + if (ib < isml) { + kb = isml - ib; + } else if (ib > ibig) { + kb = ibig - ib; + } + } + if (kb != 0) { +/* Computing MAX */ + i__2 = -kb; + i__1 = f2cmax(i__2,0); + r__1 = myscale_(&c_b29, &i__1); + i__3 = f2cmax(kb,0); + r__2 = myscale_(&c_b29, &i__3); + zlascl_("G", &c__0, &c__0, &r__1, &r__2, p, n, &b[b_offset], ldb, + info); + } + + i__1 = f2cmax(1,*m); + cnrm = zlange_("M", m, &c__1, &c__[1], &i__1, rwork); + i__1 = f2cmax(1,*p); + dnrm = zlange_("M", p, &c__1, &d__[1], &i__1, rwork); + ic = -myhuge_(&c__0); + if (cnrm > 0. && cnrm <= myhuge_(&c_b25)) { + ic = myexp_(&cnrm) + ka; + } + id = -myhuge_(&c__0); + if (dnrm > 0. && dnrm <= myhuge_(&c_b25)) { + id = myexp_(&dnrm) + kb; + } + iu = f2cmax(ic,id); + kt = 0; + if (iu > ibig) { + kt = ibig - iu; + } else if (iu < isml && (real) iu != -myhuge_(&c__0)) { + kt = isml - iu; + } + if (ka + kt != 0) { + i__1 = ka + kt; + r__1 = myscale_(&c_b29, &i__1); + zdscal_(m, &r__1, &c__[1], &c__1); + } + if (kb + kt != 0) { + i__1 = kb + kt; + r__1 = myscale_(&c_b29, &i__1); + zdscal_(p, &r__1, &d__[1], &c__1); + } + /* Compute the GRQ factorization of matrices B and A: */ /* B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P */ @@ -600,7 +715,7 @@ f"> */ /* Update c1 */ i__1 = *n - *p; - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; zgemv_("No transpose", &i__1, p, &z__1, &a[(*n - *p + 1) * a_dim1 + 1] , lda, &d__[1], &c__1, &c_b1, &c__[1], &c__1); } @@ -630,7 +745,7 @@ f"> */ nr = *m + *p - *n; if (nr > 0) { i__1 = *n - *m; - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; zgemv_("No transpose", &nr, &i__1, &z__1, &a[*n - *p + 1 + (*m + 1) * a_dim1], lda, &d__[nr + 1], &c__1, &c_b1, &c__[*n - * p + 1], &c__1); @@ -641,7 +756,7 @@ f"> */ if (nr > 0) { ztrmv_("Upper", "No transpose", "Non unit", &nr, &a[*n - *p + 1 + (*n - *p + 1) * a_dim1], lda, &d__[1], &c__1); - z__1.r = -1., z__1.i = 0.; + z__1.r = -1., z__1.i = -0.; zaxpy_(&nr, &z__1, &d__[1], &c__1, &c__[*n - *p + 1], &c__1); } @@ -650,6 +765,22 @@ f"> */ i__1 = *lwork - *p - mn; zunmrq_("Left", "Conjugate Transpose", n, &c__1, p, &b[b_offset], ldb, & work[1], &x[1], n, &work[*p + mn + 1], &i__1, info); + +/* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) */ +/* carries 2**(KA+KT) */ + + if (kt != 0) { + i__1 = -kt; + r__1 = myscale_(&c_b29, &i__1); + zdscal_(n, &r__1, &x[1], &c__1); + } + if (ka + kt != 0 && *m + *p > *n) { + i__1 = *m + *p - *n; + i__2 = -ka - kt; + r__1 = myscale_(&c_b29, &i__2); + zdscal_(&i__1, &r__1, &c__[*n - *p + 1], &c__1); + } + /* Computing MAX */ i__4 = *p + mn + 1; i__2 = lopt, i__3 = (integer) work[i__4].r; diff --git a/lapack-netlib/SRC/zgglse.f b/lapack-netlib/SRC/zgglse.f index cc558bc407..1966061ca5 100644 --- a/lapack-netlib/SRC/zgglse.f +++ b/lapack-netlib/SRC/zgglse.f @@ -5,7 +5,6 @@ * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * -*> \htmlonly *> Download ZGGLSE + dependencies *> *> [TGZ] @@ -13,7 +12,6 @@ *> [ZIP] *> *> [TXT] -*> \endhtmlonly * * Definition: * =========== @@ -52,6 +50,16 @@ *> matrices (B, A) given by *> *> B = (0 R)*Q, A = Z*T*Q. +*> +*> Callers of this subroutine should note that the singularity/rank-deficiency checks +*> implemented in this subroutine are rudimentary. The ZTRTRS subroutine called by this +*> subroutine only signals a failure due to singularity if the problem is exactly singular. +*> +*> It is conceivable for one (or more) of the factors involved in the generalized RQ +*> factorization of the pair (B, A) to be subnormally close to singularity without this +*> subroutine signalling an error. The solutions computed for such almost-rank-deficient +*> problems may be less accurate due to a loss of numerical precision. +*> *> \endverbatim * * Arguments: @@ -153,12 +161,12 @@ *> = 0: successful exit. *> < 0: if INFO = -i, the i-th argument had an illegal value. *> = 1: the upper triangular factor R associated with B in the -*> generalized RQ factorization of the pair (B, A) is +*> generalized RQ factorization of the pair (B, A) is exactly *> singular, so that rank(B) < P; the least squares *> solution could not be computed. *> = 2: the (N-P) by (N-P) part of the upper trapezoidal factor *> T associated with A in the generalized RQ factorization -*> of the pair (B, A) is singular, so that +*> of the pair (B, A) is exactly singular, so that *> rank( (A) ) < N; the least squares solution could not *> ( (B) ) *> be computed. @@ -172,11 +180,13 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \ingroup complex16OTHERsolve +*> \ingroup gglse * * ===================================================================== - SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, + SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, + $ LWORK, $ INFO ) + IMPLICIT NONE * * -- LAPACK driver routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -193,6 +203,8 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * ===================================================================== * * .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) COMPLEX*16 CONE PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) ) * .. @@ -200,17 +212,23 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, LOGICAL LQUERY INTEGER LOPT, LWKMIN, LWKOPT, MN, NB, NB1, NB2, NB3, $ NB4, NR + INTEGER IA, IB, IBIG, IC, ID, ISML, IU, KA, KB, KT + DOUBLE PRECISION ANRM, BIGNUM, BNRM, CNRM, DNRM, SMLNUM +* .. +* .. Local Arrays .. + DOUBLE PRECISION RWORK( 1 ) * .. * .. External Subroutines .. - EXTERNAL XERBLA, ZAXPY, ZCOPY, ZGEMV, ZGGRQF, ZTRMV, - $ ZTRTRS, ZUNMQR, ZUNMRQ + EXTERNAL XERBLA, ZAXPY, ZCOPY, ZDSCAL, ZGEMV, ZGGRQF, + $ ZLASCL, ZTRMV, ZTRTRS, ZUNMQR, ZUNMRQ * .. * .. External Functions .. INTEGER ILAENV - EXTERNAL ILAENV + DOUBLE PRECISION DLAMCH, ZLANGE + EXTERNAL DLAMCH, ILAENV, ZLANGE * .. * .. Intrinsic Functions .. - INTRINSIC INT, MAX, MIN + INTRINSIC EXPONENT, HUGE, INT, MAX, MIN, SCALE * .. * .. Executable Statements .. * @@ -265,6 +283,72 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( N.EQ.0 ) $ RETURN * +* Get machine parameters +* + SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' ) + BIGNUM = ONE / SMLNUM + ISML = EXPONENT( SMLNUM ) + IBIG = EXPONENT( BIGNUM ) - 1 +* +* Scale A, B, c and d by powers of two so that their largest +* entries lie in [SMLNUM,BIGNUM), the values whose EXPONENT lies +* in [ISML,IBIG]. A and B are scaled independently, by 2**KA and +* 2**KB. c is scaled by 2**(KA+KT) and d by 2**(KB+KT), which +* keeps the constraint and the residual consistent and scales x by +* 2**KT; KT brings the larger of the two right-hand sides into +* range. Scaling by a power of two is exact. A norm that is zero, +* infinite or NaN takes no part. ZLASCL is called with both +* endpoints at or above one so that forming the factor raises no +* underflow. +* + ANRM = ZLANGE( 'M', M, N, A, LDA, RWORK ) + KA = 0 + IF( ANRM.GT.ZERO .AND. ANRM.LE.HUGE( ZERO ) ) THEN + IA = EXPONENT( ANRM ) + IF( IA.LT.ISML ) THEN + KA = ISML - IA + ELSE IF( IA.GT.IBIG ) THEN + KA = IBIG - IA + END IF + END IF + IF( KA.NE.0 ) + $ CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KA, 0 ) ), + $ SCALE( ONE, MAX( KA, 0 ) ), M, N, A, LDA, INFO ) +* + BNRM = ZLANGE( 'M', P, N, B, LDB, RWORK ) + KB = 0 + IF( BNRM.GT.ZERO .AND. BNRM.LE.HUGE( ZERO ) ) THEN + IB = EXPONENT( BNRM ) + IF( IB.LT.ISML ) THEN + KB = ISML - IB + ELSE IF( IB.GT.IBIG ) THEN + KB = IBIG - IB + END IF + END IF + IF( KB.NE.0 ) + $ CALL ZLASCL( 'G', 0, 0, SCALE( ONE, MAX( -KB, 0 ) ), + $ SCALE( ONE, MAX( KB, 0 ) ), P, N, B, LDB, INFO ) +* + CNRM = ZLANGE( 'M', M, 1, C, MAX( 1, M ), RWORK ) + DNRM = ZLANGE( 'M', P, 1, D, MAX( 1, P ), RWORK ) + IC = -HUGE( 0 ) + IF( CNRM.GT.ZERO .AND. CNRM.LE.HUGE( ZERO ) ) + $ IC = EXPONENT( CNRM ) + KA + ID = -HUGE( 0 ) + IF( DNRM.GT.ZERO .AND. DNRM.LE.HUGE( ZERO ) ) + $ ID = EXPONENT( DNRM ) + KB + IU = MAX( IC, ID ) + KT = 0 + IF( IU.GT.IBIG ) THEN + KT = IBIG - IU + ELSE IF( IU.LT.ISML .AND. IU.NE.-HUGE( 0 ) ) THEN + KT = ISML - IU + END IF + IF( KA+KT.NE.0 ) + $ CALL ZDSCAL( M, SCALE( ONE, KA+KT ), C, 1 ) + IF( KB+KT.NE.0 ) + $ CALL ZDSCAL( P, SCALE( ONE, KB+KT ), D, 1 ) +* * Compute the GRQ factorization of matrices B and A: * * B*Q**H = ( 0 T12 ) P Z**H*A*Q**H = ( R11 R12 ) N-P @@ -303,7 +387,8 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * * Update c1 * - CALL ZGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), LDA, + CALL ZGEMV( 'No transpose', N-P, P, -CONE, A( 1, N-P+1 ), + $ LDA, $ D, 1, CONE, C, 1 ) END IF * @@ -328,7 +413,8 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, IF( M.LT.N ) THEN NR = M + P - N IF( NR.GT.0 ) - $ CALL ZGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, M+1 ), + $ CALL ZGEMV( 'No transpose', NR, N-M, -CONE, A( N-P+1, + $ M+1 ), $ LDA, D( NR+1 ), 1, CONE, C( N-P+1 ), 1 ) ELSE NR = P @@ -343,6 +429,15 @@ SUBROUTINE ZGGLSE( M, N, P, A, LDA, B, LDB, C, D, X, WORK, LWORK, * CALL ZUNMRQ( 'Left', 'Conjugate Transpose', N, 1, P, B, LDB, $ WORK( 1 ), X, N, WORK( P+MN+1 ), LWORK-P-MN, INFO ) +* +* Undo scaling: x carries 2**KT and the residual in c(N-P+1:M) +* carries 2**(KA+KT) +* + IF( KT.NE.0 ) + $ CALL ZDSCAL( N, SCALE( ONE, -KT ), X, 1 ) + IF( KA+KT.NE.0 .AND. M+P.GT.N ) + $ CALL ZDSCAL( M+P-N, SCALE( ONE, -KA-KT ), C( N-P+1 ), 1 ) +* WORK( 1 ) = P + MN + MAX( LOPT, INT( WORK( P+MN+1 ) ) ) * RETURN diff --git a/lapack-netlib/TESTING/EIG/cglmts.f b/lapack-netlib/TESTING/EIG/cglmts.f index ecf8699f27..8445bcb32c 100644 --- a/lapack-netlib/TESTING/EIG/cglmts.f +++ b/lapack-netlib/TESTING/EIG/cglmts.f @@ -28,7 +28,7 @@ *> \verbatim *> *> CGLMTS tests CGGGLM - a subroutine for solving the generalized -*> linear model problem. +*> linear model problem, including independent scaling of its inputs. *> \endverbatim * * Arguments: @@ -147,6 +147,7 @@ * ===================================================================== SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, $ X, U, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -167,22 +168,28 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, * .. Parameters .. REAL ZERO PARAMETER ( ZERO = 0.0E+0 ) + REAL ONE + PARAMETER ( ONE = 1.0E+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) COMPLEX CONE PARAMETER ( CONE = 1.0E+0 ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL + REAL ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. REAL SCASUM, SLAMCH, CLANGE - EXTERNAL SCASUM, SLAMCH, CLANGE + LOGICAL SISNAN + EXTERNAL SISNAN, SCASUM, SLAMCH, CLANGE * .. * .. External Subroutines .. - EXTERNAL CLACPY + EXTERNAL CCOPY, CGEMV, CGGGLM, CLACPY, CSSCAL * * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -225,6 +232,68 @@ SUBROUTINE CGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF +* +* A, B and d may be scaled independently: for factors a, b, d, +* the solutions become (d/a)*x and (d/b)*u. Test large and +* tiny inputs, then undo solution scaling before the residual. +* + TNRM = MAX( CLANGE( 'M', N, M, A, LDA, RWORK ), + $ CLANGE( 'M', N, P, B, LDB, RWORK ), + $ CLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny d, tiny (A,d), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = SLAMCH( 'Safe minimum' ) / + $ SLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + CALL CLACPY( 'Full', N, M, A, LDA, AF, LDA ) + CALL CLACPY( 'Full', N, P, B, LDB, BF, LDB ) + CALL CCOPY( N, D, 1, DF, 1 ) + DO 10 J = 1, M + CALL CSSCAL( N, ASCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL CSSCAL( N, BSCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL CSSCAL( N, DSCL, DF, 1 ) +* + CALL CGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT = ONE / EPS + GO TO 30 + END IF + CALL CSSCAL( M, ASCL / DSCL, X, 1 ) + CALL CSSCAL( P, BSCL / DSCL, U, 1 ) +* + CALL CCOPY( N, D, 1, DF, 1 ) + CALL CGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, + $ DF, 1 ) + CALL CGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, + $ DF, 1 ) + DNORM = SCASUM( N, DF, 1 ) + XNORM = SCASUM( M, X, 1 ) + SCASUM( P, U, 1 ) + IF( SISNAN( DNORM ) .OR. SISNAN( XNORM ) ) THEN + RESULT = ONE / EPS + ELSE IF( XNORM.GT.ZERO ) THEN + RESULT = MAX( RESULT, + $ ( ( DNORM / YNORM ) / XNORM ) / EPS ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/clsets.f b/lapack-netlib/TESTING/EIG/clsets.f index 2124a9e091..fad4a0169f 100644 --- a/lapack-netlib/TESTING/EIG/clsets.f +++ b/lapack-netlib/TESTING/EIG/clsets.f @@ -27,7 +27,7 @@ *> \verbatim *> *> CLSETS tests CGGLSE - a subroutine for solving linear equality -*> constrained least square problem (LSE). +*> constrained least square problem (LSE), including scaled inputs. *> \endverbatim * * Arguments: @@ -152,6 +152,7 @@ * ===================================================================== SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, $ D, DF, X, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -170,10 +171,25 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J + REAL ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM +* .. +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) +* .. +* .. External Functions .. + LOGICAL SISNAN + REAL CLANGE, SLAMCH + EXTERNAL SISNAN, CLANGE, SLAMCH * .. * .. External Subroutines .. - EXTERNAL CGGLSE, CLACPY, CGET02 + EXTERNAL CCOPY, CGGLSE, CLACPY, CGET02, CSSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -203,6 +219,76 @@ SUBROUTINE CLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * CALL CGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, $ RWORK, RESULT( 2 ) ) +* +* Scaling (A,c) and (B,d) independently leaves x unchanged. +* Scaling only (c,d) scales x by the same factor. Exercise +* both ends of the range and check residuals at the input scale. +* + TNRM = MAX( CLANGE( 'M', M, N, A, LDA, RWORK ), + $ CLANGE( 'M', P, N, B, LDB, RWORK ), + $ CLANGE( 'M', M, 1, C, M, RWORK ), + $ CLANGE( 'M', P, 1, D, P, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny (c,d), tiny (A,c), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = SLAMCH( 'Safe minimum' ) / + $ SLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + CSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + IF( ISCALE.EQ.4 ) DSCL = ONE + IF( ISCALE.EQ.5 ) CSCL = ONE + CALL CLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL CLACPY( 'Full', P, N, B, LDB, BF, LDB ) + CALL CCOPY( M, C, 1, CF, 1 ) + CALL CCOPY( P, D, 1, DF, 1 ) + DO 10 J = 1, N + CALL CSSCAL( M, ASCL, AF( 1, J ), 1 ) + CALL CSSCAL( P, BSCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL CSSCAL( M, CSCL, CF, 1 ) + CALL CSSCAL( P, DSCL, DF, 1 ) +* + CALL CGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + RESULT( 2 ) = RESULT( 1 ) + GO TO 30 + END IF + IF( ISCALE.EQ.3 ) + $ CALL CSSCAL( N, ONE / SCL, X, 1 ) +* + CALL CCOPY( M, C, 1, CF, 1 ) + CALL CCOPY( P, D, 1, DF, 1 ) + CALL CGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M, + $ RWORK, RESID ) + IF( SISNAN( RESID ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = MAX( RESULT( 1 ), RESID ) + END IF + CALL CGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, + $ RWORK, RESID ) + IF( SISNAN( RESID ) ) THEN + RESULT( 2 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 2 ) = MAX( RESULT( 2 ), RESID ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/dglmts.f b/lapack-netlib/TESTING/EIG/dglmts.f index 0286a757cd..03be4917a1 100644 --- a/lapack-netlib/TESTING/EIG/dglmts.f +++ b/lapack-netlib/TESTING/EIG/dglmts.f @@ -24,7 +24,7 @@ *> \verbatim *> *> DGLMTS tests DGGGLM - a subroutine for solving the generalized -*> linear model problem. +*> linear model problem, including independent scaling of its inputs. *> \endverbatim * * Arguments: @@ -143,6 +143,7 @@ * ===================================================================== SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, $ WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -163,21 +164,25 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM + DOUBLE PRECISION ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. DOUBLE PRECISION DASUM, DLAMCH, DLANGE - EXTERNAL DASUM, DLAMCH, DLANGE + LOGICAL DISNAN + EXTERNAL DISNAN, DASUM, DLAMCH, DLANGE * .. * .. External Subroutines .. * - EXTERNAL DCOPY, DGEMV, DGGGLM, DLACPY + EXTERNAL DCOPY, DGEMV, DGGGLM, DLACPY, DSCAL * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -218,6 +223,68 @@ SUBROUTINE DGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF +* +* A, B and d may be scaled independently: for factors a, b, d, +* the solutions become (d/a)*x and (d/b)*u. Test large and +* tiny inputs, then undo solution scaling before the residual. +* + TNRM = MAX( DLANGE( 'M', N, M, A, LDA, RWORK ), + $ DLANGE( 'M', N, P, B, LDB, RWORK ), + $ DLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny d, tiny (A,d), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = DLAMCH( 'Safe minimum' ) / + $ DLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + CALL DLACPY( 'Full', N, M, A, LDA, AF, LDA ) + CALL DLACPY( 'Full', N, P, B, LDB, BF, LDB ) + CALL DCOPY( N, D, 1, DF, 1 ) + DO 10 J = 1, M + CALL DSCAL( N, ASCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL DSCAL( N, BSCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL DSCAL( N, DSCL, DF, 1 ) +* + CALL DGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT = ONE / EPS + GO TO 30 + END IF + CALL DSCAL( M, ASCL / DSCL, X, 1 ) + CALL DSCAL( P, BSCL / DSCL, U, 1 ) +* + CALL DCOPY( N, D, 1, DF, 1 ) + CALL DGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, + $ DF, 1 ) + CALL DGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, + $ DF, 1 ) + DNORM = DASUM( N, DF, 1 ) + XNORM = DASUM( M, X, 1 ) + DASUM( P, U, 1 ) + IF( DISNAN( DNORM ) .OR. DISNAN( XNORM ) ) THEN + RESULT = ONE / EPS + ELSE IF( XNORM.GT.ZERO ) THEN + RESULT = MAX( RESULT, + $ ( ( DNORM / YNORM ) / XNORM ) / EPS ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/dlsets.f b/lapack-netlib/TESTING/EIG/dlsets.f index 3616183a52..4ec6926b84 100644 --- a/lapack-netlib/TESTING/EIG/dlsets.f +++ b/lapack-netlib/TESTING/EIG/dlsets.f @@ -23,7 +23,7 @@ *> \verbatim *> *> DLSETS tests DGGLSE - a subroutine for solving linear equality -*> constrained least square problem (LSE). +*> constrained least square problem (LSE), including scaled inputs. *> \endverbatim * * Arguments: @@ -148,6 +148,7 @@ * ===================================================================== SUBROUTINE DLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, $ X, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -165,10 +166,25 @@ SUBROUTINE DLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, $ RESULT( 2 ), RWORK( * ), WORK( LWORK ), X( * ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J + DOUBLE PRECISION ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM +* .. +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) +* .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION DLANGE, DLAMCH + EXTERNAL DISNAN, DLANGE, DLAMCH * .. * .. External Subroutines .. - EXTERNAL DCOPY, DGET02, DGGLSE, DLACPY + EXTERNAL DCOPY, DGET02, DGGLSE, DLACPY, DSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -198,6 +214,76 @@ SUBROUTINE DLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, * CALL DGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, RWORK, $ RESULT( 2 ) ) +* +* Scaling (A,c) and (B,d) independently leaves x unchanged. +* Scaling only (c,d) scales x by the same factor. Exercise +* both ends of the range and check residuals at the input scale. +* + TNRM = MAX( DLANGE( 'M', M, N, A, LDA, RWORK ), + $ DLANGE( 'M', P, N, B, LDB, RWORK ), + $ DLANGE( 'M', M, 1, C, M, RWORK ), + $ DLANGE( 'M', P, 1, D, P, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny (c,d), tiny (A,c), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = DLAMCH( 'Safe minimum' ) / + $ DLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + CSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + IF( ISCALE.EQ.4 ) DSCL = ONE + IF( ISCALE.EQ.5 ) CSCL = ONE + CALL DLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL DLACPY( 'Full', P, N, B, LDB, BF, LDB ) + CALL DCOPY( M, C, 1, CF, 1 ) + CALL DCOPY( P, D, 1, DF, 1 ) + DO 10 J = 1, N + CALL DSCAL( M, ASCL, AF( 1, J ), 1 ) + CALL DSCAL( P, BSCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL DSCAL( M, CSCL, CF, 1 ) + CALL DSCAL( P, DSCL, DF, 1 ) +* + CALL DGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + RESULT( 2 ) = RESULT( 1 ) + GO TO 30 + END IF + IF( ISCALE.EQ.3 ) + $ CALL DSCAL( N, ONE / SCL, X, 1 ) +* + CALL DCOPY( M, C, 1, CF, 1 ) + CALL DCOPY( P, D, 1, DF, 1 ) + CALL DGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M, + $ RWORK, RESID ) + IF( DISNAN( RESID ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = MAX( RESULT( 1 ), RESID ) + END IF + CALL DGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, + $ RWORK, RESID ) + IF( DISNAN( RESID ) ) THEN + RESULT( 2 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 2 ) = MAX( RESULT( 2 ), RESID ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/sglmts.f b/lapack-netlib/TESTING/EIG/sglmts.f index 0859b3e0dc..1766745300 100644 --- a/lapack-netlib/TESTING/EIG/sglmts.f +++ b/lapack-netlib/TESTING/EIG/sglmts.f @@ -27,7 +27,7 @@ *> \verbatim *> *> SGLMTS tests SGGGLM - a subroutine for solving the generalized -*> linear model problem. +*> linear model problem, including independent scaling of its inputs. *> \endverbatim * * Arguments: @@ -146,6 +146,7 @@ * ===================================================================== SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, $ X, U, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -165,20 +166,24 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, * .. Parameters .. REAL ZERO, ONE PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J REAL ANORM, BNORM, EPS, XNORM, YNORM, DNORM, UNFL + REAL ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. REAL SASUM, SLAMCH, SLANGE - EXTERNAL SASUM, SLAMCH, SLANGE + LOGICAL SISNAN + EXTERNAL SISNAN, SASUM, SLAMCH, SLANGE * .. * .. External Subroutines .. - EXTERNAL SLACPY + EXTERNAL SCOPY, SGEMV, SGGGLM, SLACPY, SSCAL * * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -221,6 +226,68 @@ SUBROUTINE SGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) /EPS END IF +* +* A, B and d may be scaled independently: for factors a, b, d, +* the solutions become (d/a)*x and (d/b)*u. Test large and +* tiny inputs, then undo solution scaling before the residual. +* + TNRM = MAX( SLANGE( 'M', N, M, A, LDA, RWORK ), + $ SLANGE( 'M', N, P, B, LDB, RWORK ), + $ SLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny d, tiny (A,d), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = SLAMCH( 'Safe minimum' ) / + $ SLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + CALL SLACPY( 'Full', N, M, A, LDA, AF, LDA ) + CALL SLACPY( 'Full', N, P, B, LDB, BF, LDB ) + CALL SCOPY( N, D, 1, DF, 1 ) + DO 10 J = 1, M + CALL SSCAL( N, ASCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL SSCAL( N, BSCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL SSCAL( N, DSCL, DF, 1 ) +* + CALL SGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT = ONE / EPS + GO TO 30 + END IF + CALL SSCAL( M, ASCL / DSCL, X, 1 ) + CALL SSCAL( P, BSCL / DSCL, U, 1 ) +* + CALL SCOPY( N, D, 1, DF, 1 ) + CALL SGEMV( 'No transpose', N, M, -ONE, A, LDA, X, 1, ONE, + $ DF, 1 ) + CALL SGEMV( 'No transpose', N, P, -ONE, B, LDB, U, 1, ONE, + $ DF, 1 ) + DNORM = SASUM( N, DF, 1 ) + XNORM = SASUM( M, X, 1 ) + SASUM( P, U, 1 ) + IF( SISNAN( DNORM ) .OR. SISNAN( XNORM ) ) THEN + RESULT = ONE / EPS + ELSE IF( XNORM.GT.ZERO ) THEN + RESULT = MAX( RESULT, + $ ( ( DNORM / YNORM ) / XNORM ) / EPS ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/slsets.f b/lapack-netlib/TESTING/EIG/slsets.f index 0fb66decf2..6dffa75048 100644 --- a/lapack-netlib/TESTING/EIG/slsets.f +++ b/lapack-netlib/TESTING/EIG/slsets.f @@ -27,7 +27,7 @@ *> \verbatim *> *> SLSETS tests SGGLSE - a subroutine for solving linear equality -*> constrained least square problem (LSE). +*> constrained least square problem (LSE), including scaled inputs. *> \endverbatim * * Arguments: @@ -152,6 +152,7 @@ * ===================================================================== SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, $ D, DF, X, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -170,10 +171,25 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J + REAL ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM +* .. +* .. Parameters .. + REAL ZERO, ONE + PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) +* .. +* .. External Functions .. + LOGICAL SISNAN + REAL SLANGE, SLAMCH + EXTERNAL SISNAN, SLANGE, SLAMCH * .. * .. External Subroutines .. - EXTERNAL SGGLSE, SLACPY, SGET02 + EXTERNAL SCOPY, SGGLSE, SLACPY, SGET02, SSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -203,6 +219,76 @@ SUBROUTINE SLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, * CALL SGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, $ RWORK, RESULT( 2 ) ) +* +* Scaling (A,c) and (B,d) independently leaves x unchanged. +* Scaling only (c,d) scales x by the same factor. Exercise +* both ends of the range and check residuals at the input scale. +* + TNRM = MAX( SLANGE( 'M', M, N, A, LDA, RWORK ), + $ SLANGE( 'M', P, N, B, LDB, RWORK ), + $ SLANGE( 'M', M, 1, C, M, RWORK ), + $ SLANGE( 'M', P, 1, D, P, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.SLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny (c,d), tiny (A,c), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = SLAMCH( 'Safe minimum' ) / + $ SLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + CSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + IF( ISCALE.EQ.4 ) DSCL = ONE + IF( ISCALE.EQ.5 ) CSCL = ONE + CALL SLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL SLACPY( 'Full', P, N, B, LDB, BF, LDB ) + CALL SCOPY( M, C, 1, CF, 1 ) + CALL SCOPY( P, D, 1, DF, 1 ) + DO 10 J = 1, N + CALL SSCAL( M, ASCL, AF( 1, J ), 1 ) + CALL SSCAL( P, BSCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL SSCAL( M, CSCL, CF, 1 ) + CALL SSCAL( P, DSCL, DF, 1 ) +* + CALL SGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + RESULT( 2 ) = RESULT( 1 ) + GO TO 30 + END IF + IF( ISCALE.EQ.3 ) + $ CALL SSCAL( N, ONE / SCL, X, 1 ) +* + CALL SCOPY( M, C, 1, CF, 1 ) + CALL SCOPY( P, D, 1, DF, 1 ) + CALL SGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M, + $ RWORK, RESID ) + IF( SISNAN( RESID ) ) THEN + RESULT( 1 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = MAX( RESULT( 1 ), RESID ) + END IF + CALL SGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, + $ RWORK, RESID ) + IF( SISNAN( RESID ) ) THEN + RESULT( 2 ) = ONE / SLAMCH( 'Epsilon' ) + ELSE + RESULT( 2 ) = MAX( RESULT( 2 ), RESID ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/zglmts.f b/lapack-netlib/TESTING/EIG/zglmts.f index 63664241d3..7649808693 100644 --- a/lapack-netlib/TESTING/EIG/zglmts.f +++ b/lapack-netlib/TESTING/EIG/zglmts.f @@ -24,7 +24,7 @@ *> \verbatim *> *> ZGLMTS tests ZGGGLM - a subroutine for solving the generalized -*> linear model problem. +*> linear model problem, including independent scaling of its inputs. *> \endverbatim * * Arguments: @@ -143,6 +143,7 @@ * ===================================================================== SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, $ WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -164,23 +165,29 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, * .. Parameters .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D+0 ) + DOUBLE PRECISION ONE + PARAMETER ( ONE = 1.0D+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) COMPLEX*16 CONE PARAMETER ( CONE = 1.0D+0 ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J DOUBLE PRECISION ANORM, BNORM, DNORM, EPS, UNFL, XNORM, YNORM + DOUBLE PRECISION ASCL, BSCL, DSCL, SCL, TNRM * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DZASUM, ZLANGE - EXTERNAL DLAMCH, DZASUM, ZLANGE + LOGICAL DISNAN + EXTERNAL DISNAN, DLAMCH, DZASUM, ZLANGE * .. * .. External Subroutines .. * - EXTERNAL ZCOPY, ZGEMV, ZGGGLM, ZLACPY + EXTERNAL ZCOPY, ZGEMV, ZGGGLM, ZLACPY, ZDSCAL * .. * .. Intrinsic Functions .. - INTRINSIC MAX + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -223,6 +230,68 @@ SUBROUTINE ZGLMTS( N, M, P, A, AF, LDA, B, BF, LDB, D, DF, X, U, ELSE RESULT = ( ( DNORM / YNORM ) / XNORM ) / EPS END IF +* +* A, B and d may be scaled independently: for factors a, b, d, +* the solutions become (d/a)*x and (d/b)*u. Test large and +* tiny inputs, then undo solution scaling before the residual. +* + TNRM = MAX( ZLANGE( 'M', N, M, A, LDA, RWORK ), + $ ZLANGE( 'M', N, P, B, LDB, RWORK ), + $ ZLANGE( 'M', N, 1, D, N, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny d, tiny (A,d), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = DLAMCH( 'Safe minimum' ) / + $ DLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + CALL ZLACPY( 'Full', N, M, A, LDA, AF, LDA ) + CALL ZLACPY( 'Full', N, P, B, LDB, BF, LDB ) + CALL ZCOPY( N, D, 1, DF, 1 ) + DO 10 J = 1, M + CALL ZDSCAL( N, ASCL, AF( 1, J ), 1 ) + 10 CONTINUE + DO 20 J = 1, P + CALL ZDSCAL( N, BSCL, BF( 1, J ), 1 ) + 20 CONTINUE + CALL ZDSCAL( N, DSCL, DF, 1 ) +* + CALL ZGGGLM( N, M, P, AF, LDA, BF, LDB, DF, X, U, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT = ONE / EPS + GO TO 30 + END IF + CALL ZDSCAL( M, ASCL / DSCL, X, 1 ) + CALL ZDSCAL( P, BSCL / DSCL, U, 1 ) +* + CALL ZCOPY( N, D, 1, DF, 1 ) + CALL ZGEMV( 'No transpose', N, M, -CONE, A, LDA, X, 1, CONE, + $ DF, 1 ) + CALL ZGEMV( 'No transpose', N, P, -CONE, B, LDB, U, 1, CONE, + $ DF, 1 ) + DNORM = DZASUM( N, DF, 1 ) + XNORM = DZASUM( M, X, 1 ) + DZASUM( P, U, 1 ) + IF( DISNAN( DNORM ) .OR. DISNAN( XNORM ) ) THEN + RESULT = ONE / EPS + ELSE IF( XNORM.GT.ZERO ) THEN + RESULT = MAX( RESULT, + $ ( ( DNORM / YNORM ) / XNORM ) / EPS ) + END IF + 30 CONTINUE + END IF * RETURN * diff --git a/lapack-netlib/TESTING/EIG/zlsets.f b/lapack-netlib/TESTING/EIG/zlsets.f index 8c8ff01b22..c9d00ab4e1 100644 --- a/lapack-netlib/TESTING/EIG/zlsets.f +++ b/lapack-netlib/TESTING/EIG/zlsets.f @@ -23,7 +23,7 @@ *> \verbatim *> *> ZLSETS tests ZGGLSE - a subroutine for solving linear equality -*> constrained least square problem (LSE). +*> constrained least square problem (LSE), including scaled inputs. *> \endverbatim * * Arguments: @@ -148,6 +148,7 @@ * ===================================================================== SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, $ X, WORK, LWORK, RWORK, RESULT ) + IMPLICIT NONE * * -- LAPACK test routine -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- @@ -166,10 +167,25 @@ SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, $ WORK( LWORK ), X( * ) * .. * .. Local Scalars .. - INTEGER INFO + INTEGER INFO, ISCALE, J + DOUBLE PRECISION ASCL, BSCL, CSCL, DSCL, RESID, SCL, TNRM +* .. +* .. Parameters .. + DOUBLE PRECISION ZERO, ONE + PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) + INTEGER MAXEXP + PARAMETER ( MAXEXP = MAXEXPONENT( ZERO ) - 2 ) +* .. +* .. External Functions .. + LOGICAL DISNAN + DOUBLE PRECISION ZLANGE, DLAMCH + EXTERNAL DISNAN, ZLANGE, DLAMCH * .. * .. External Subroutines .. - EXTERNAL ZCOPY, ZGET02, ZGGLSE, ZLACPY + EXTERNAL ZCOPY, ZGET02, ZGGLSE, ZLACPY, ZDSCAL +* .. +* .. Intrinsic Functions .. + INTRINSIC EXPONENT, MAX, MAXEXPONENT, SCALE * .. * .. Executable Statements .. * @@ -199,6 +215,76 @@ SUBROUTINE ZLSETS( M, P, N, A, AF, LDA, B, BF, LDB, C, CF, D, DF, * CALL ZGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, RWORK, $ RESULT( 2 ) ) +* +* Scaling (A,c) and (B,d) independently leaves x unchanged. +* Scaling only (c,d) scales x by the same factor. Exercise +* both ends of the range and check residuals at the input scale. +* + TNRM = MAX( ZLANGE( 'M', M, N, A, LDA, RWORK ), + $ ZLANGE( 'M', P, N, B, LDB, RWORK ), + $ ZLANGE( 'M', M, 1, C, M, RWORK ), + $ ZLANGE( 'M', P, 1, D, P, RWORK ) ) + IF( TNRM.GT.ZERO .AND. TNRM.LE.DLAMCH( 'Overflow' ) ) THEN +* +* Cases: common large, common tiny, tiny (c,d), tiny (A,c), +* and tiny (B,d). +* + DO 30 ISCALE = 1, 5 + IF( ISCALE.EQ.1 ) THEN + SCL = SCALE( ONE, MAXEXP-EXPONENT( TNRM ) ) + ELSE + SCL = DLAMCH( 'Safe minimum' ) / + $ DLAMCH( 'Precision' ) + SCL = SCALE( ONE, EXPONENT( SCL )-4- + $ EXPONENT( TNRM ) ) + END IF + ASCL = SCL + BSCL = SCL + CSCL = SCL + DSCL = SCL + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.5 ) ASCL = ONE + IF( ISCALE.EQ.3 .OR. ISCALE.EQ.4 ) BSCL = ONE + IF( ISCALE.EQ.4 ) DSCL = ONE + IF( ISCALE.EQ.5 ) CSCL = ONE + CALL ZLACPY( 'Full', M, N, A, LDA, AF, LDA ) + CALL ZLACPY( 'Full', P, N, B, LDB, BF, LDB ) + CALL ZCOPY( M, C, 1, CF, 1 ) + CALL ZCOPY( P, D, 1, DF, 1 ) + DO 10 J = 1, N + CALL ZDSCAL( M, ASCL, AF( 1, J ), 1 ) + CALL ZDSCAL( P, BSCL, BF( 1, J ), 1 ) + 10 CONTINUE + CALL ZDSCAL( M, CSCL, CF, 1 ) + CALL ZDSCAL( P, DSCL, DF, 1 ) +* + CALL ZGGLSE( M, N, P, AF, LDA, BF, LDB, CF, DF, X, WORK, + $ LWORK, INFO ) + IF( INFO.NE.0 ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + RESULT( 2 ) = RESULT( 1 ) + GO TO 30 + END IF + IF( ISCALE.EQ.3 ) + $ CALL ZDSCAL( N, ONE / SCL, X, 1 ) +* + CALL ZCOPY( M, C, 1, CF, 1 ) + CALL ZCOPY( P, D, 1, DF, 1 ) + CALL ZGET02( 'No transpose', M, N, 1, A, LDA, X, N, CF, M, + $ RWORK, RESID ) + IF( DISNAN( RESID ) ) THEN + RESULT( 1 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 1 ) = MAX( RESULT( 1 ), RESID ) + END IF + CALL ZGET02( 'No transpose', P, N, 1, B, LDB, X, N, DF, P, + $ RWORK, RESID ) + IF( DISNAN( RESID ) ) THEN + RESULT( 2 ) = ONE / DLAMCH( 'Epsilon' ) + ELSE + RESULT( 2 ) = MAX( RESULT( 2 ), RESID ) + END IF + 30 CONTINUE + END IF * RETURN *