diff --git a/ext/f2c_blas/idamax.c b/ext/f2c_blas/idamax.c deleted file mode 100644 index 3037a85bc..000000000 --- a/ext/f2c_blas/idamax.c +++ /dev/null @@ -1,67 +0,0 @@ -#include "blaswrap.h" -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -integer idamax_(integer *n, doublereal *dx, integer *incx) -{ - /* System generated locals */ - integer ret_val, i__1; - doublereal d__1; - /* Local variables */ - static doublereal dmax__; - static integer i__, ix; -/* finds the index of element having max. absolute value. - jack dongarra, linpack, 3/11/78. - modified 3/93 to return if incx .le. 0. - modified 12/3/93, array(1) declarations changed to array(*) - Parameter adjustments */ - --dx; - /* Function Body */ - ret_val = 0; - if (*n < 1 || *incx <= 0) { - return ret_val; - } - ret_val = 1; - if (*n == 1) { - return ret_val; - } - if (*incx == 1) { - goto L20; - } -/* code for increment not equal to 1 */ - ix = 1; - dmax__ = abs(dx[1]); - ix += *incx; - i__1 = *n; - for (i__ = 2; i__ <= i__1; ++i__) { - if ((d__1 = dx[ix], abs(d__1)) <= dmax__) { - goto L5; - } - ret_val = i__; - dmax__ = (d__1 = dx[ix], abs(d__1)); -L5: - ix += *incx; -/* L10: */ - } - return ret_val; -/* code for increment equal to 1 */ -L20: - dmax__ = abs(dx[1]); - i__1 = *n; - for (i__ = 2; i__ <= i__1; ++i__) { - if ((d__1 = dx[i__], abs(d__1)) <= dmax__) { - goto L30; - } - ret_val = i__; - dmax__ = (d__1 = dx[i__], abs(d__1)); -L30: - ; - } - return ret_val; -} /* idamax_ */ - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_lapack/dgbsvx.c b/ext/f2c_lapack/dgbsvx.c deleted file mode 100644 index 36ff0add1..000000000 --- a/ext/f2c_lapack/dgbsvx.c +++ /dev/null @@ -1,642 +0,0 @@ -#include "blaswrap.h" -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -/* Subroutine */ int dgbsvx_(char *fact, char *trans, integer *n, integer *kl, - integer *ku, integer *nrhs, doublereal *ab, integer *ldab, - doublereal *afb, integer *ldafb, integer *ipiv, char *equed, - doublereal *r__, doublereal *c__, doublereal *b, integer *ldb, - doublereal *x, integer *ldx, doublereal *rcond, doublereal *ferr, - doublereal *berr, doublereal *work, integer *iwork, integer *info) -{ -/* -- LAPACK driver routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - June 30, 1999 - - - Purpose - ======= - - DGBSVX uses the LU factorization to compute the solution to a real - system of linear equations A * X = B, A**T * X = B, or A**H * X = B, - where A is a band matrix of order N with KL subdiagonals and KU - superdiagonals, and X and B are N-by-NRHS matrices. - - Error bounds on the solution and a condition estimate are also - provided. - - Description - =========== - - The following steps are performed by this subroutine: - - 1. If FACT = 'E', real scaling factors are computed to equilibrate - the system: - TRANS = 'N': diag(R)*A*diag(C) *inv(diag(C))*X = diag(R)*B - TRANS = 'T': (diag(R)*A*diag(C))**T *inv(diag(R))*X = diag(C)*B - TRANS = 'C': (diag(R)*A*diag(C))**H *inv(diag(R))*X = diag(C)*B - Whether or not the system will be equilibrated depends on the - scaling of the matrix A, but if equilibration is used, A is - overwritten by diag(R)*A*diag(C) and B by diag(R)*B (if TRANS='N') - or diag(C)*B (if TRANS = 'T' or 'C'). - - 2. If FACT = 'N' or 'E', the LU decomposition is used to factor the - matrix A (after equilibration if FACT = 'E') as - A = L * U, - where L is a product of permutation and unit lower triangular - matrices with KL subdiagonals, and U is upper triangular with - KL+KU superdiagonals. - - 3. If some U(i,i)=0, so that U is exactly singular, then the routine - returns with INFO = i. Otherwise, the factored form of A is used - to estimate the condition number of the matrix A. If the - reciprocal of the condition number is less than machine precision, - INFO = N+1 is returned as a warning, but the routine still goes on - to solve for X and compute error bounds as described below. - - 4. The system of equations is solved for X using the factored form - of A. - - 5. Iterative refinement is applied to improve the computed solution - matrix and calculate error bounds and backward error estimates - for it. - - 6. If equilibration was used, the matrix X is premultiplied by - diag(C) (if TRANS = 'N') or diag(R) (if TRANS = 'T' or 'C') so - that it solves the original system before equilibration. - - Arguments - ========= - - FACT (input) CHARACTER*1 - Specifies whether or not the factored form of the matrix A is - supplied on entry, and if not, whether the matrix A should be - equilibrated before it is factored. - = 'F': On entry, AFB and IPIV contain the factored form of - A. If EQUED is not 'N', the matrix A has been - equilibrated with scaling factors given by R and C. - AB, AFB, and IPIV are not modified. - = 'N': The matrix A will be copied to AFB and factored. - = 'E': The matrix A will be equilibrated if necessary, then - copied to AFB and factored. - - TRANS (input) CHARACTER*1 - Specifies the form of the system of equations. - = 'N': A * X = B (No transpose) - = 'T': A**T * X = B (Transpose) - = 'C': A**H * X = B (Transpose) - - N (input) INTEGER - The number of linear equations, i.e., the order of the - matrix A. N >= 0. - - KL (input) INTEGER - The number of subdiagonals within the band of A. KL >= 0. - - KU (input) INTEGER - The number of superdiagonals within the band of A. KU >= 0. - - NRHS (input) INTEGER - The number of right hand sides, i.e., the number of columns - of the matrices B and X. NRHS >= 0. - - AB (input/output) DOUBLE PRECISION array, dimension (LDAB,N) - On entry, the matrix A in band storage, in rows 1 to KL+KU+1. - The j-th column of A is stored in the j-th column of the - array AB as follows: - AB(KU+1+i-j,j) = A(i,j) for max(1,j-KU)<=i<=min(N,j+kl) - - If FACT = 'F' and EQUED is not 'N', then A must have been - equilibrated by the scaling factors in R and/or C. AB is not - modified if FACT = 'F' or 'N', or if FACT = 'E' and - EQUED = 'N' on exit. - - On exit, if EQUED .ne. 'N', A is scaled as follows: - EQUED = 'R': A := diag(R) * A - EQUED = 'C': A := A * diag(C) - EQUED = 'B': A := diag(R) * A * diag(C). - - LDAB (input) INTEGER - The leading dimension of the array AB. LDAB >= KL+KU+1. - - AFB (input or output) DOUBLE PRECISION array, dimension (LDAFB,N) - If FACT = 'F', then AFB is an input argument and on entry - contains details of the LU factorization of the band matrix - A, as computed by DGBTRF. U is stored as an upper triangular - band matrix with KL+KU superdiagonals in rows 1 to KL+KU+1, - and the multipliers used during the factorization are stored - in rows KL+KU+2 to 2*KL+KU+1. If EQUED .ne. 'N', then AFB is - the factored form of the equilibrated matrix A. - - If FACT = 'N', then AFB is an output argument and on exit - returns details of the LU factorization of A. - - If FACT = 'E', then AFB is an output argument and on exit - returns details of the LU factorization of the equilibrated - matrix A (see the description of AB for the form of the - equilibrated matrix). - - LDAFB (input) INTEGER - The leading dimension of the array AFB. LDAFB >= 2*KL+KU+1. - - IPIV (input or output) INTEGER array, dimension (N) - If FACT = 'F', then IPIV is an input argument and on entry - contains the pivot indices from the factorization A = L*U - as computed by DGBTRF; row i of the matrix was interchanged - with row IPIV(i). - - If FACT = 'N', then IPIV is an output argument and on exit - contains the pivot indices from the factorization A = L*U - of the original matrix A. - - If FACT = 'E', then IPIV is an output argument and on exit - contains the pivot indices from the factorization A = L*U - of the equilibrated matrix A. - - EQUED (input or output) CHARACTER*1 - Specifies the form of equilibration that was done. - = 'N': No equilibration (always true if FACT = 'N'). - = 'R': Row equilibration, i.e., A has been premultiplied by - diag(R). - = 'C': Column equilibration, i.e., A has been postmultiplied - by diag(C). - = 'B': Both row and column equilibration, i.e., A has been - replaced by diag(R) * A * diag(C). - EQUED is an input argument if FACT = 'F'; otherwise, it is an - output argument. - - R (input or output) DOUBLE PRECISION array, dimension (N) - The row scale factors for A. If EQUED = 'R' or 'B', A is - multiplied on the left by diag(R); if EQUED = 'N' or 'C', R - is not accessed. R is an input argument if FACT = 'F'; - otherwise, R is an output argument. If FACT = 'F' and - EQUED = 'R' or 'B', each element of R must be positive. - - C (input or output) DOUBLE PRECISION array, dimension (N) - The column scale factors for A. If EQUED = 'C' or 'B', A is - multiplied on the right by diag(C); if EQUED = 'N' or 'R', C - is not accessed. C is an input argument if FACT = 'F'; - otherwise, C is an output argument. If FACT = 'F' and - EQUED = 'C' or 'B', each element of C must be positive. - - B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) - On entry, the right hand side matrix B. - On exit, - if EQUED = 'N', B is not modified; - if TRANS = 'N' and EQUED = 'R' or 'B', B is overwritten by - diag(R)*B; - if TRANS = 'T' or 'C' and EQUED = 'C' or 'B', B is - overwritten by diag(C)*B. - - LDB (input) INTEGER - The leading dimension of the array B. LDB >= max(1,N). - - X (output) DOUBLE PRECISION array, dimension (LDX,NRHS) - If INFO = 0 or INFO = N+1, the N-by-NRHS solution matrix X - to the original system of equations. Note that A and B are - modified on exit if EQUED .ne. 'N', and the solution to the - equilibrated system is inv(diag(C))*X if TRANS = 'N' and - EQUED = 'C' or 'B', or inv(diag(R))*X if TRANS = 'T' or 'C' - and EQUED = 'R' or 'B'. - - LDX (input) INTEGER - The leading dimension of the array X. LDX >= max(1,N). - - RCOND (output) DOUBLE PRECISION - The estimate of the reciprocal condition number of the matrix - A after equilibration (if done). If RCOND is less than the - machine precision (in particular, if RCOND = 0), the matrix - is singular to working precision. This condition is - indicated by a return code of INFO > 0. - - FERR (output) DOUBLE PRECISION array, dimension (NRHS) - The estimated forward error bound for each solution vector - X(j) (the j-th column of the solution matrix X). - If XTRUE is the true solution corresponding to X(j), FERR(j) - is an estimated upper bound for the magnitude of the largest - element in (X(j) - XTRUE) divided by the magnitude of the - largest element in X(j). The estimate is as reliable as - the estimate for RCOND, and is almost always a slight - overestimate of the true error. - - BERR (output) DOUBLE PRECISION array, dimension (NRHS) - The componentwise relative backward error of each solution - vector X(j) (i.e., the smallest relative change in - any element of A or B that makes X(j) an exact solution). - - WORK (workspace/output) DOUBLE PRECISION array, dimension (3*N) - On exit, WORK(1) contains the reciprocal pivot growth - factor norm(A)/norm(U). The "max absolute element" norm is - used. If WORK(1) is much less than 1, then the stability - of the LU factorization of the (equilibrated) matrix A - could be poor. This also means that the solution X, condition - estimator RCOND, and forward error bound FERR could be - unreliable. If factorization fails with 0 0: if INFO = i, and i is - <= N: U(i,i) is exactly zero. The factorization - has been completed, but the factor U is exactly - singular, so the solution and error bounds - could not be computed. RCOND = 0 is returned. - = N+1: U is nonsingular, but RCOND is less than machine - precision, meaning that the matrix is singular - to working precision. Nevertheless, the - solution and error bounds are computed because - there are a number of situations where the - computed solution can be more accurate than the - value of RCOND would suggest. - - ===================================================================== - - - Parameter adjustments */ - /* Table of constant values */ - static integer c__1 = 1; - - /* System generated locals */ - integer ab_dim1, ab_offset, afb_dim1, afb_offset, b_dim1, b_offset, - x_dim1, x_offset, i__1, i__2, i__3, i__4, i__5; - doublereal d__1, d__2, d__3; - /* Local variables */ - static doublereal amax; - static char norm[1]; - static integer i__, j; - extern logical lsame_(char *, char *); - static doublereal rcmin, rcmax, anorm; - extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, - doublereal *, integer *); - static logical equil; - static integer j1, j2; - extern doublereal dlamch_(char *), dlangb_(char *, integer *, - integer *, integer *, doublereal *, integer *, doublereal *); - extern /* Subroutine */ int dlaqgb_(integer *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - doublereal *, doublereal *, doublereal *, char *), - dgbcon_(char *, integer *, integer *, integer *, doublereal *, - integer *, integer *, doublereal *, doublereal *, doublereal *, - integer *, integer *); - static doublereal colcnd; - extern doublereal dlantb_(char *, char *, char *, integer *, integer *, - doublereal *, integer *, doublereal *); - extern /* Subroutine */ int dgbequ_(integer *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - doublereal *, doublereal *, doublereal *, integer *), dgbrfs_( - char *, integer *, integer *, integer *, integer *, doublereal *, - integer *, doublereal *, integer *, integer *, doublereal *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - doublereal *, integer *, integer *), dgbtrf_(integer *, - integer *, integer *, integer *, doublereal *, integer *, integer - *, integer *); - static logical nofact; - extern /* Subroutine */ int dlacpy_(char *, integer *, integer *, - doublereal *, integer *, doublereal *, integer *), - xerbla_(char *, integer *); - static doublereal bignum; - extern /* Subroutine */ int dgbtrs_(char *, integer *, integer *, integer - *, integer *, doublereal *, integer *, integer *, doublereal *, - integer *, integer *); - static integer infequ; - static logical colequ; - static doublereal rowcnd; - static logical notran; - static doublereal smlnum; - static logical rowequ; - static doublereal rpvgrw; -#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] -#define x_ref(a_1,a_2) x[(a_2)*x_dim1 + a_1] -#define ab_ref(a_1,a_2) ab[(a_2)*ab_dim1 + a_1] -#define afb_ref(a_1,a_2) afb[(a_2)*afb_dim1 + a_1] - - - ab_dim1 = *ldab; - ab_offset = 1 + ab_dim1 * 1; - ab -= ab_offset; - afb_dim1 = *ldafb; - afb_offset = 1 + afb_dim1 * 1; - afb -= afb_offset; - --ipiv; - --r__; - --c__; - b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; - b -= b_offset; - x_dim1 = *ldx; - x_offset = 1 + x_dim1 * 1; - x -= x_offset; - --ferr; - --berr; - --work; - --iwork; - - /* Function Body */ - *info = 0; - nofact = lsame_(fact, "N"); - equil = lsame_(fact, "E"); - notran = lsame_(trans, "N"); - if (nofact || equil) { - *(unsigned char *)equed = 'N'; - rowequ = FALSE_; - colequ = FALSE_; - } else { - rowequ = lsame_(equed, "R") || lsame_(equed, - "B"); - colequ = lsame_(equed, "C") || lsame_(equed, - "B"); - smlnum = dlamch_("Safe minimum"); - bignum = 1. / smlnum; - } - -/* Test the input parameters. */ - - if (! nofact && ! equil && ! lsame_(fact, "F")) { - *info = -1; - } else if (! notran && ! lsame_(trans, "T") && ! - lsame_(trans, "C")) { - *info = -2; - } else if (*n < 0) { - *info = -3; - } else if (*kl < 0) { - *info = -4; - } else if (*ku < 0) { - *info = -5; - } else if (*nrhs < 0) { - *info = -6; - } else if (*ldab < *kl + *ku + 1) { - *info = -8; - } else if (*ldafb < (*kl << 1) + *ku + 1) { - *info = -10; - } else if (lsame_(fact, "F") && ! (rowequ || colequ - || lsame_(equed, "N"))) { - *info = -12; - } else { - if (rowequ) { - rcmin = bignum; - rcmax = 0.; - i__1 = *n; - for (j = 1; j <= i__1; ++j) { -/* Computing MIN */ - d__1 = rcmin, d__2 = r__[j]; - rcmin = min(d__1,d__2); -/* Computing MAX */ - d__1 = rcmax, d__2 = r__[j]; - rcmax = max(d__1,d__2); -/* L10: */ - } - if (rcmin <= 0.) { - *info = -13; - } else if (*n > 0) { - rowcnd = max(rcmin,smlnum) / min(rcmax,bignum); - } else { - rowcnd = 1.; - } - } - if (colequ && *info == 0) { - rcmin = bignum; - rcmax = 0.; - i__1 = *n; - for (j = 1; j <= i__1; ++j) { -/* Computing MIN */ - d__1 = rcmin, d__2 = c__[j]; - rcmin = min(d__1,d__2); -/* Computing MAX */ - d__1 = rcmax, d__2 = c__[j]; - rcmax = max(d__1,d__2); -/* L20: */ - } - if (rcmin <= 0.) { - *info = -14; - } else if (*n > 0) { - colcnd = max(rcmin,smlnum) / min(rcmax,bignum); - } else { - colcnd = 1.; - } - } - if (*info == 0) { - if (*ldb < max(1,*n)) { - *info = -16; - } else if (*ldx < max(1,*n)) { - *info = -18; - } - } - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("DGBSVX", &i__1); - return 0; - } - - if (equil) { - -/* Compute row and column scalings to equilibrate the matrix A. */ - - dgbequ_(n, n, kl, ku, &ab[ab_offset], ldab, &r__[1], &c__[1], &rowcnd, - &colcnd, &amax, &infequ); - if (infequ == 0) { - -/* Equilibrate the matrix. */ - - dlaqgb_(n, n, kl, ku, &ab[ab_offset], ldab, &r__[1], &c__[1], & - rowcnd, &colcnd, &amax, equed); - rowequ = lsame_(equed, "R") || lsame_(equed, - "B"); - colequ = lsame_(equed, "C") || lsame_(equed, - "B"); - } - } - -/* Scale the right hand side. */ - - if (notran) { - if (rowequ) { - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - b_ref(i__, j) = r__[i__] * b_ref(i__, j); -/* L30: */ - } -/* L40: */ - } - } - } else if (colequ) { - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - b_ref(i__, j) = c__[i__] * b_ref(i__, j); -/* L50: */ - } -/* L60: */ - } - } - - if (nofact || equil) { - -/* Compute the LU factorization of the band matrix A. */ - - i__1 = *n; - for (j = 1; j <= i__1; ++j) { -/* Computing MAX */ - i__2 = j - *ku; - j1 = max(i__2,1); -/* Computing MIN */ - i__2 = j + *kl; - j2 = min(i__2,*n); - i__2 = j2 - j1 + 1; - dcopy_(&i__2, &ab_ref(*ku + 1 - j + j1, j), &c__1, &afb_ref(*kl + - *ku + 1 - j + j1, j), &c__1); -/* L70: */ - } - - dgbtrf_(n, n, kl, ku, &afb[afb_offset], ldafb, &ipiv[1], info); - -/* Return if INFO is non-zero. */ - - if (*info != 0) { - if (*info > 0) { - -/* Compute the reciprocal pivot growth factor of the - leading rank-deficient INFO columns of A. */ - - anorm = 0.; - i__1 = *info; - for (j = 1; j <= i__1; ++j) { -/* Computing MAX */ - i__2 = *ku + 2 - j; -/* Computing MIN */ - i__4 = *n + *ku + 1 - j, i__5 = *kl + *ku + 1; - i__3 = min(i__4,i__5); - for (i__ = max(i__2,1); i__ <= i__3; ++i__) { -/* Computing MAX */ - d__2 = anorm, d__3 = (d__1 = ab_ref(i__, j), abs(d__1) - ); - anorm = max(d__2,d__3); -/* L80: */ - } -/* L90: */ - } -/* Computing MAX */ - i__1 = 1, i__3 = *kl + *ku + 2 - *info; -/* Computing MIN */ - i__4 = *info - 1, i__5 = *kl + *ku; - i__2 = min(i__4,i__5); - rpvgrw = dlantb_("M", "U", "N", info, &i__2, &afb_ref(max( - i__1,i__3), 1), ldafb, &work[1]); - if (rpvgrw == 0.) { - rpvgrw = 1.; - } else { - rpvgrw = anorm / rpvgrw; - } - work[1] = rpvgrw; - *rcond = 0.; - } - return 0; - } - } - -/* Compute the norm of the matrix A and the - reciprocal pivot growth factor RPVGRW. */ - - if (notran) { - *(unsigned char *)norm = '1'; - } else { - *(unsigned char *)norm = 'I'; - } - anorm = dlangb_(norm, n, kl, ku, &ab[ab_offset], ldab, &work[1]); - i__1 = *kl + *ku; - rpvgrw = dlantb_("M", "U", "N", n, &i__1, &afb[afb_offset], ldafb, &work[ - 1]); - if (rpvgrw == 0.) { - rpvgrw = 1.; - } else { - rpvgrw = dlangb_("M", n, kl, ku, &ab[ab_offset], ldab, &work[1]) / rpvgrw; - } - -/* Compute the reciprocal of the condition number of A. */ - - dgbcon_(norm, n, kl, ku, &afb[afb_offset], ldafb, &ipiv[1], &anorm, rcond, - &work[1], &iwork[1], info); - -/* Set INFO = N+1 if the matrix is singular to working precision. */ - - if (*rcond < dlamch_("Epsilon")) { - *info = *n + 1; - } - -/* Compute the solution matrix X. */ - - dlacpy_("Full", n, nrhs, &b[b_offset], ldb, &x[x_offset], ldx); - dgbtrs_(trans, n, kl, ku, nrhs, &afb[afb_offset], ldafb, &ipiv[1], &x[ - x_offset], ldx, info); - -/* Use iterative refinement to improve the computed solution and - compute error bounds and backward error estimates for it. */ - - dgbrfs_(trans, n, kl, ku, nrhs, &ab[ab_offset], ldab, &afb[afb_offset], - ldafb, &ipiv[1], &b[b_offset], ldb, &x[x_offset], ldx, &ferr[1], & - berr[1], &work[1], &iwork[1], info); - -/* Transform the solution matrix X to a solution of the original - system. */ - - if (notran) { - if (colequ) { - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__3 = *n; - for (i__ = 1; i__ <= i__3; ++i__) { - x_ref(i__, j) = c__[i__] * x_ref(i__, j); -/* L100: */ - } -/* L110: */ - } - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - ferr[j] /= colcnd; -/* L120: */ - } - } - } else if (rowequ) { - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__3 = *n; - for (i__ = 1; i__ <= i__3; ++i__) { - x_ref(i__, j) = r__[i__] * x_ref(i__, j); -/* L130: */ - } -/* L140: */ - } - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - ferr[j] /= rowcnd; -/* L150: */ - } - } - - work[1] = rpvgrw; - return 0; - -/* End of DGBSVX */ - -} /* dgbsvx_ */ - -#undef afb_ref -#undef ab_ref -#undef x_ref -#undef b_ref - - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_lapack/dgelsd.c b/ext/f2c_lapack/dgelsd.c deleted file mode 100644 index 3718ca05e..000000000 --- a/ext/f2c_lapack/dgelsd.c +++ /dev/null @@ -1,684 +0,0 @@ -#include "blaswrap.h" -/* -- translated by f2c (version 19990503). - You must link the resulting object file with the libraries: - -lf2c -lm (in that order) -*/ - -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -/* Table of constant values */ - -static integer c__6 = 6; -static integer c_n1 = -1; -static integer c__9 = 9; -static integer c__0 = 0; -static integer c__1 = 1; -static doublereal c_b82 = 0.; - -/* Subroutine */ int dgelsd_(integer *m, integer *n, integer *nrhs, - doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * - s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, - integer *iwork, integer *info) -{ - /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2, i__3, i__4; - - /* Builtin functions */ - double log(doublereal); - - /* Local variables */ - static doublereal anrm, bnrm; - static integer itau, nlvl, iascl, ibscl; - static doublereal sfmin; - static integer minmn, maxmn, itaup, itauq, mnthr, nwork; - extern /* Subroutine */ int dlabad_(doublereal *, doublereal *); - static integer ie, il; - extern /* Subroutine */ int dgebrd_(integer *, integer *, doublereal *, - integer *, doublereal *, doublereal *, doublereal *, doublereal *, - doublereal *, integer *, integer *); - extern doublereal dlamch_(char *); - static integer mm; - extern doublereal dlange_(char *, integer *, integer *, doublereal *, - integer *, doublereal *); - extern /* Subroutine */ int dgelqf_(integer *, integer *, doublereal *, - integer *, doublereal *, doublereal *, integer *, integer *), - dlalsd_(char *, integer *, integer *, integer *, doublereal *, - doublereal *, doublereal *, integer *, doublereal *, integer *, - doublereal *, integer *, integer *), dlascl_(char *, - integer *, integer *, doublereal *, doublereal *, integer *, - integer *, doublereal *, integer *, integer *), dgeqrf_( - integer *, integer *, doublereal *, integer *, doublereal *, - doublereal *, integer *, integer *), dlacpy_(char *, integer *, - integer *, doublereal *, integer *, doublereal *, integer *), dlaset_(char *, integer *, integer *, doublereal *, - doublereal *, doublereal *, integer *), xerbla_(char *, - integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - static doublereal bignum; - extern /* Subroutine */ int dormbr_(char *, char *, char *, integer *, - integer *, integer *, doublereal *, integer *, doublereal *, - doublereal *, integer *, doublereal *, integer *, integer *); - static integer wlalsd; - extern /* Subroutine */ int dormlq_(char *, char *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - integer *, doublereal *, integer *, integer *); - static integer ldwork; - extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - integer *, doublereal *, integer *, integer *); - static integer minwrk, maxwrk; - static doublereal smlnum; - static logical lquery; - static integer smlsiz; - static doublereal eps; - - -#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] -#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] - - -/* -- LAPACK driver routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - October 31, 1999 - - - Purpose - ======= - - DGELSD computes the minimum-norm solution to a real linear least - squares problem: - minimize 2-norm(| b - A*x |) - using the singular value decomposition (SVD) of A. A is an M-by-N - matrix which may be rank-deficient. - - Several right hand side vectors b and solution vectors x can be - handled in a single call; they are stored as the columns of the - M-by-NRHS right hand side matrix B and the N-by-NRHS solution - matrix X. - - The problem is solved in three steps: - (1) Reduce the coefficient matrix A to bidiagonal form with - Householder transformations, reducing the original problem - into a "bidiagonal least squares problem" (BLS) - (2) Solve the BLS using a divide and conquer approach. - (3) Apply back all the Householder tranformations to solve - the original least squares problem. - - The effective rank of A is determined by treating as zero those - singular values which are less than RCOND times the largest singular - value. - - The divide and conquer algorithm makes very mild assumptions about - floating point arithmetic. It will work on machines with a guard - digit in add/subtract, or on those binary machines without guard - digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or - Cray-2. It could conceivably fail on hexadecimal or decimal machines - without guard digits, but we know of none. - - Arguments - ========= - - M (input) INTEGER - The number of rows of A. M >= 0. - - N (input) INTEGER - The number of columns of A. N >= 0. - - NRHS (input) INTEGER - The number of right hand sides, i.e., the number of columns - of the matrices B and X. NRHS >= 0. - - A (input) DOUBLE PRECISION array, dimension (LDA,N) - On entry, the M-by-N matrix A. - On exit, A has been destroyed. - - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,M). - - B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) - On entry, the M-by-NRHS right hand side matrix B. - On exit, B is overwritten by the N-by-NRHS solution - matrix X. If m >= n and RANK = n, the residual - sum-of-squares for the solution in the i-th column is given - by the sum of squares of elements n+1:m in that column. - - LDB (input) INTEGER - The leading dimension of the array B. LDB >= max(1,max(M,N)). - - S (output) DOUBLE PRECISION array, dimension (min(M,N)) - The singular values of A in decreasing order. - The condition number of A in the 2-norm = S(1)/S(min(m,n)). - - RCOND (input) DOUBLE PRECISION - RCOND is used to determine the effective rank of A. - Singular values S(i) <= RCOND*S(1) are treated as zero. - If RCOND < 0, machine precision is used instead. - - RANK (output) INTEGER - The effective rank of A, i.e., the number of singular values - which are greater than RCOND*S(1). - - WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) - On exit, if INFO = 0, WORK(1) returns the optimal LWORK. - - LWORK (input) INTEGER - The dimension of the array WORK. LWORK must be at least 1. - The exact minimum amount of workspace needed depends on M, - N and NRHS. As long as LWORK is at least - 12*N + 2*N*SMLSIZ + 8*N*NLVL + N*NRHS + (SMLSIZ+1)**2, - if M is greater than or equal to N or - 12*M + 2*M*SMLSIZ + 8*M*NLVL + M*NRHS + (SMLSIZ+1)**2, - if M is less than N, the code will execute correctly. - SMLSIZ is returned by ILAENV and is equal to the maximum - size of the subproblems at the bottom of the computation - tree (usually about 25), and - NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 ) - For good performance, LWORK should generally be larger. - - If LWORK = -1, then a workspace query is assumed; the routine - only calculates the optimal size of the WORK array, returns - this value as the first entry of the WORK array, and no error - message related to LWORK is issued by XERBLA. - - IWORK (workspace) INTEGER array, dimension (LIWORK) - LIWORK >= 3 * MINMN * NLVL + 11 * MINMN, - where MINMN = MIN( M,N ). - - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value. - > 0: the algorithm for computing the SVD failed to converge; - if INFO = i, i off-diagonal elements of an intermediate - bidiagonal form did not converge to zero. - - Further Details - =============== - - Based on contributions by - Ming Gu and Ren-Cang Li, Computer Science Division, University of - California at Berkeley, USA - Osni Marques, LBNL/NERSC, USA - - ===================================================================== - - - Test the input arguments. - - Parameter adjustments */ - a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; - a -= a_offset; - b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; - b -= b_offset; - --s; - --work; - --iwork; - - /* Function Body */ - *info = 0; - minmn = min(*m,*n); - maxmn = max(*m,*n); - mnthr = ilaenv_(&c__6, "DGELSD", " ", m, n, nrhs, &c_n1, (ftnlen)6, ( - ftnlen)1); - lquery = *lwork == -1; - if (*m < 0) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*nrhs < 0) { - *info = -3; - } else if (*lda < max(1,*m)) { - *info = -5; - } else if (*ldb < max(1,maxmn)) { - *info = -7; - } - - smlsiz = ilaenv_(&c__9, "DGELSD", " ", &c__0, &c__0, &c__0, &c__0, ( - ftnlen)6, (ftnlen)1); - -/* Compute workspace. - (Note: Comments in the code beginning "Workspace:" describe the - minimal amount of workspace needed at that point in the code, - as well as the preferred amount for good performance. - NB refers to the optimal block size for the immediately - following subroutine, as returned by ILAENV.) */ - - minwrk = 1; - minmn = max(1,minmn); -/* Computing MAX */ - i__1 = (integer) (log((doublereal) minmn / (doublereal) (smlsiz + 1)) / - log(2.)) + 1; - nlvl = max(i__1,0); - - if (*info == 0) { - maxwrk = 0; - mm = *m; - if (*m >= *n && *m >= mnthr) { - -/* Path 1a - overdetermined, with many more rows than columns. */ - - mm = *n; -/* Computing MAX */ - i__1 = maxwrk, i__2 = *n + *n * ilaenv_(&c__1, "DGEQRF", " ", m, - n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *n + *nrhs * ilaenv_(&c__1, "DORMQR", "LT", - m, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)2); - maxwrk = max(i__1,i__2); - } - if (*m >= *n) { - -/* Path 1 - overdetermined or exactly determined. - - Computing MAX */ - i__1 = maxwrk, i__2 = *n * 3 + (mm + *n) * ilaenv_(&c__1, "DGEBRD" - , " ", &mm, n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *n * 3 + *nrhs * ilaenv_(&c__1, "DORMBR", - "QLT", &mm, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *n * 3 + (*n - 1) * ilaenv_(&c__1, "DORMBR", - "PLN", n, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); - maxwrk = max(i__1,i__2); -/* Computing 2nd power */ - i__1 = smlsiz + 1; - wlalsd = *n * 9 + (*n << 1) * smlsiz + (*n << 3) * nlvl + *n * * - nrhs + i__1 * i__1; -/* Computing MAX */ - i__1 = maxwrk, i__2 = *n * 3 + wlalsd; - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = *n * 3 + mm, i__2 = *n * 3 + *nrhs, i__1 = max(i__1,i__2), - i__2 = *n * 3 + wlalsd; - minwrk = max(i__1,i__2); - } - if (*n > *m) { -/* Computing 2nd power */ - i__1 = smlsiz + 1; - wlalsd = *m * 9 + (*m << 1) * smlsiz + (*m << 3) * nlvl + *m * * - nrhs + i__1 * i__1; - if (*n >= mnthr) { - -/* Path 2a - underdetermined, with many more columns - than rows. */ - - maxwrk = *m + *m * ilaenv_(&c__1, "DGELQF", " ", m, n, &c_n1, - &c_n1, (ftnlen)6, (ftnlen)1); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m << 1) * - ilaenv_(&c__1, "DGEBRD", " ", m, m, &c_n1, &c_n1, ( - ftnlen)6, (ftnlen)1); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + *nrhs * ilaenv_(& - c__1, "DORMBR", "QLT", m, nrhs, m, &c_n1, (ftnlen)6, ( - ftnlen)3); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + (*m - 1) * - ilaenv_(&c__1, "DORMBR", "PLN", m, nrhs, m, &c_n1, ( - ftnlen)6, (ftnlen)3); - maxwrk = max(i__1,i__2); - if (*nrhs > 1) { -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + *m + *m * *nrhs; - maxwrk = max(i__1,i__2); - } else { -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + (*m << 1); - maxwrk = max(i__1,i__2); - } -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m + *nrhs * ilaenv_(&c__1, "DORMLQ", - "LT", n, nrhs, m, &c_n1, (ftnlen)6, (ftnlen)2); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * *m + (*m << 2) + wlalsd; - maxwrk = max(i__1,i__2); - } else { - -/* Path 2 - remaining underdetermined cases. */ - - maxwrk = *m * 3 + (*n + *m) * ilaenv_(&c__1, "DGEBRD", " ", m, - n, &c_n1, &c_n1, (ftnlen)6, (ftnlen)1); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * 3 + *nrhs * ilaenv_(&c__1, "DORMBR" - , "QLT", m, nrhs, n, &c_n1, (ftnlen)6, (ftnlen)3); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * 3 + *m * ilaenv_(&c__1, "DORMBR", - "PLN", n, nrhs, m, &c_n1, (ftnlen)6, (ftnlen)3); - maxwrk = max(i__1,i__2); -/* Computing MAX */ - i__1 = maxwrk, i__2 = *m * 3 + wlalsd; - maxwrk = max(i__1,i__2); - } -/* Computing MAX */ - i__1 = *m * 3 + *nrhs, i__2 = *m * 3 + *m, i__1 = max(i__1,i__2), - i__2 = *m * 3 + wlalsd; - minwrk = max(i__1,i__2); - } - minwrk = min(minwrk,maxwrk); - work[1] = (doublereal) maxwrk; - if (*lwork < minwrk && ! lquery) { - *info = -12; - } - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("DGELSD", &i__1); - return 0; - } else if (lquery) { - goto L10; - } - -/* Quick return if possible. */ - - if (*m == 0 || *n == 0) { - *rank = 0; - return 0; - } - -/* Get machine parameters. */ - - eps = dlamch_("P"); - sfmin = dlamch_("S"); - smlnum = sfmin / eps; - bignum = 1. / smlnum; - dlabad_(&smlnum, &bignum); - -/* Scale A if max entry outside range [SMLNUM,BIGNUM]. */ - - anrm = dlange_("M", m, n, &a[a_offset], lda, &work[1]); - iascl = 0; - if (anrm > 0. && anrm < smlnum) { - -/* Scale matrix norm up to SMLNUM. */ - - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, - info); - iascl = 1; - } else if (anrm > bignum) { - -/* Scale matrix norm down to BIGNUM. */ - - dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, - info); - iascl = 2; - } else if (anrm == 0.) { - -/* Matrix all zero. Return zero solution. */ - - i__1 = max(*m,*n); - dlaset_("F", &i__1, nrhs, &c_b82, &c_b82, &b[b_offset], ldb); - dlaset_("F", &minmn, &c__1, &c_b82, &c_b82, &s[1], &c__1); - *rank = 0; - goto L10; - } - -/* Scale B if max entry outside range [SMLNUM,BIGNUM]. */ - - bnrm = dlange_("M", m, nrhs, &b[b_offset], ldb, &work[1]); - ibscl = 0; - if (bnrm > 0. && bnrm < smlnum) { - -/* Scale matrix norm up to SMLNUM. */ - - dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 1; - } else if (bnrm > bignum) { - -/* Scale matrix norm down to BIGNUM. */ - - dlascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 2; - } - -/* If M < N make sure certain entries of B are zero. */ - - if (*m < *n) { - i__1 = *n - *m; - dlaset_("F", &i__1, nrhs, &c_b82, &c_b82, &b_ref(*m + 1, 1), ldb); - } - -/* Overdetermined case. */ - - if (*m >= *n) { - -/* Path 1 - overdetermined or exactly determined. */ - - mm = *m; - if (*m >= mnthr) { - -/* Path 1a - overdetermined, with many more rows than columns. */ - - mm = *n; - itau = 1; - nwork = itau + *n; - -/* Compute A=Q*R. - (Workspace: need 2*N, prefer N+N*NB) */ - - i__1 = *lwork - nwork + 1; - dgeqrf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &i__1, - info); - -/* Multiply B by transpose(Q). - (Workspace: need N+NRHS, prefer N+NRHS*NB) */ - - i__1 = *lwork - nwork + 1; - dormqr_("L", "T", m, nrhs, n, &a[a_offset], lda, &work[itau], &b[ - b_offset], ldb, &work[nwork], &i__1, info); - -/* Zero out below R. */ - - if (*n > 1) { - i__1 = *n - 1; - i__2 = *n - 1; - dlaset_("L", &i__1, &i__2, &c_b82, &c_b82, &a_ref(2, 1), lda); - } - } - - ie = 1; - itauq = ie + *n; - itaup = itauq + *n; - nwork = itaup + *n; - -/* Bidiagonalize R in A. - (Workspace: need 3*N+MM, prefer 3*N+(MM+N)*NB) */ - - i__1 = *lwork - nwork + 1; - dgebrd_(&mm, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], & - work[itaup], &work[nwork], &i__1, info); - -/* Multiply B by transpose of left bidiagonalizing vectors of R. - (Workspace: need 3*N+NRHS, prefer 3*N+NRHS*NB) */ - - i__1 = *lwork - nwork + 1; - dormbr_("Q", "L", "T", &mm, nrhs, n, &a[a_offset], lda, &work[itauq], - &b[b_offset], ldb, &work[nwork], &i__1, info); - -/* Solve the bidiagonal least squares problem. */ - - dlalsd_("U", &smlsiz, n, nrhs, &s[1], &work[ie], &b[b_offset], ldb, - rcond, rank, &work[nwork], &iwork[1], info); - if (*info != 0) { - goto L10; - } - -/* Multiply B by right bidiagonalizing vectors of R. */ - - i__1 = *lwork - nwork + 1; - dormbr_("P", "L", "N", n, nrhs, n, &a[a_offset], lda, &work[itaup], & - b[b_offset], ldb, &work[nwork], &i__1, info); - - } else /* if(complicated condition) */ { -/* Computing MAX */ - i__1 = *m, i__2 = (*m << 1) - 4, i__1 = max(i__1,i__2), i__1 = max( - i__1,*nrhs), i__2 = *n - *m * 3; - if (*n >= mnthr && *lwork >= (*m << 2) + *m * *m + max(i__1,i__2)) { - -/* Path 2a - underdetermined, with many more columns than rows - and sufficient workspace for an efficient algorithm. */ - - ldwork = *m; -/* Computing MAX - Computing MAX */ - i__3 = *m, i__4 = (*m << 1) - 4, i__3 = max(i__3,i__4), i__3 = - max(i__3,*nrhs), i__4 = *n - *m * 3; - i__1 = (*m << 2) + *m * *lda + max(i__3,i__4), i__2 = *m * *lda + - *m + *m * *nrhs; - if (*lwork >= max(i__1,i__2)) { - ldwork = *lda; - } - itau = 1; - nwork = *m + 1; - -/* Compute A=L*Q. - (Workspace: need 2*M, prefer M+M*NB) */ - - i__1 = *lwork - nwork + 1; - dgelqf_(m, n, &a[a_offset], lda, &work[itau], &work[nwork], &i__1, - info); - il = nwork; - -/* Copy L to WORK(IL), zeroing out above its diagonal. */ - - dlacpy_("L", m, m, &a[a_offset], lda, &work[il], &ldwork); - i__1 = *m - 1; - i__2 = *m - 1; - dlaset_("U", &i__1, &i__2, &c_b82, &c_b82, &work[il + ldwork], & - ldwork); - ie = il + ldwork * *m; - itauq = ie + *m; - itaup = itauq + *m; - nwork = itaup + *m; - -/* Bidiagonalize L in WORK(IL). - (Workspace: need M*M+5*M, prefer M*M+4*M+2*M*NB) */ - - i__1 = *lwork - nwork + 1; - dgebrd_(m, m, &work[il], &ldwork, &s[1], &work[ie], &work[itauq], - &work[itaup], &work[nwork], &i__1, info); - -/* Multiply B by transpose of left bidiagonalizing vectors of L. - (Workspace: need M*M+4*M+NRHS, prefer M*M+4*M+NRHS*NB) */ - - i__1 = *lwork - nwork + 1; - dormbr_("Q", "L", "T", m, nrhs, m, &work[il], &ldwork, &work[ - itauq], &b[b_offset], ldb, &work[nwork], &i__1, info); - -/* Solve the bidiagonal least squares problem. */ - - dlalsd_("U", &smlsiz, m, nrhs, &s[1], &work[ie], &b[b_offset], - ldb, rcond, rank, &work[nwork], &iwork[1], info); - if (*info != 0) { - goto L10; - } - -/* Multiply B by right bidiagonalizing vectors of L. */ - - i__1 = *lwork - nwork + 1; - dormbr_("P", "L", "N", m, nrhs, m, &work[il], &ldwork, &work[ - itaup], &b[b_offset], ldb, &work[nwork], &i__1, info); - -/* Zero out below first M rows of B. */ - - i__1 = *n - *m; - dlaset_("F", &i__1, nrhs, &c_b82, &c_b82, &b_ref(*m + 1, 1), ldb); - nwork = itau + *m; - -/* Multiply transpose(Q) by B. - (Workspace: need M+NRHS, prefer M+NRHS*NB) */ - - i__1 = *lwork - nwork + 1; - dormlq_("L", "T", n, nrhs, m, &a[a_offset], lda, &work[itau], &b[ - b_offset], ldb, &work[nwork], &i__1, info); - - } else { - -/* Path 2 - remaining underdetermined cases. */ - - ie = 1; - itauq = ie + *m; - itaup = itauq + *m; - nwork = itaup + *m; - -/* Bidiagonalize A. - (Workspace: need 3*M+N, prefer 3*M+(M+N)*NB) */ - - i__1 = *lwork - nwork + 1; - dgebrd_(m, n, &a[a_offset], lda, &s[1], &work[ie], &work[itauq], & - work[itaup], &work[nwork], &i__1, info); - -/* Multiply B by transpose of left bidiagonalizing vectors. - (Workspace: need 3*M+NRHS, prefer 3*M+NRHS*NB) */ - - i__1 = *lwork - nwork + 1; - dormbr_("Q", "L", "T", m, nrhs, n, &a[a_offset], lda, &work[itauq] - , &b[b_offset], ldb, &work[nwork], &i__1, info); - -/* Solve the bidiagonal least squares problem. */ - - dlalsd_("L", &smlsiz, m, nrhs, &s[1], &work[ie], &b[b_offset], - ldb, rcond, rank, &work[nwork], &iwork[1], info); - if (*info != 0) { - goto L10; - } - -/* Multiply B by right bidiagonalizing vectors of A. */ - - i__1 = *lwork - nwork + 1; - dormbr_("P", "L", "N", n, nrhs, m, &a[a_offset], lda, &work[itaup] - , &b[b_offset], ldb, &work[nwork], &i__1, info); - - } - } - -/* Undo scaling. */ - - if (iascl == 1) { - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("G", &c__0, &c__0, &smlnum, &anrm, &minmn, &c__1, &s[1], & - minmn, info); - } else if (iascl == 2) { - dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("G", &c__0, &c__0, &bignum, &anrm, &minmn, &c__1, &s[1], & - minmn, info); - } - if (ibscl == 1) { - dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } else if (ibscl == 2) { - dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } - -L10: - work[1] = (doublereal) maxwrk; - return 0; - -/* End of DGELSD */ - -} /* dgelsd_ */ - -#undef b_ref -#undef a_ref - - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_lapack/dgelsx.c b/ext/f2c_lapack/dgelsx.c deleted file mode 100644 index d3ac89475..000000000 --- a/ext/f2c_lapack/dgelsx.c +++ /dev/null @@ -1,424 +0,0 @@ -#include "blaswrap.h" -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -/* Subroutine */ int dgelsx_(integer *m, integer *n, integer *nrhs, - doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * - jpvt, doublereal *rcond, integer *rank, doublereal *work, integer * - info) -{ -/* -- LAPACK driver routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - March 31, 1993 - - - Purpose - ======= - - This routine is deprecated and has been replaced by routine DGELSY. - - DGELSX computes the minimum-norm solution to a real linear least - squares problem: - minimize || A * X - B || - using a complete orthogonal factorization of A. A is an M-by-N - matrix which may be rank-deficient. - - Several right hand side vectors b and solution vectors x can be - handled in a single call; they are stored as the columns of the - M-by-NRHS right hand side matrix B and the N-by-NRHS solution - matrix X. - - The routine first computes a QR factorization with column pivoting: - A * P = Q * [ R11 R12 ] - [ 0 R22 ] - with R11 defined as the largest leading submatrix whose estimated - condition number is less than 1/RCOND. The order of R11, RANK, - is the effective rank of A. - - Then, R22 is considered to be negligible, and R12 is annihilated - by orthogonal transformations from the right, arriving at the - complete orthogonal factorization: - A * P = Q * [ T11 0 ] * Z - [ 0 0 ] - The minimum-norm solution is then - X = P * Z' [ inv(T11)*Q1'*B ] - [ 0 ] - where Q1 consists of the first RANK columns of Q. - - Arguments - ========= - - M (input) INTEGER - The number of rows of the matrix A. M >= 0. - - N (input) INTEGER - The number of columns of the matrix A. N >= 0. - - NRHS (input) INTEGER - The number of right hand sides, i.e., the number of - columns of matrices B and X. NRHS >= 0. - - A (input/output) DOUBLE PRECISION array, dimension (LDA,N) - On entry, the M-by-N matrix A. - On exit, A has been overwritten by details of its - complete orthogonal factorization. - - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,M). - - B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) - On entry, the M-by-NRHS right hand side matrix B. - On exit, the N-by-NRHS solution matrix X. - If m >= n and RANK = n, the residual sum-of-squares for - the solution in the i-th column is given by the sum of - squares of elements N+1:M in that column. - - LDB (input) INTEGER - The leading dimension of the array B. LDB >= max(1,M,N). - - JPVT (input/output) INTEGER array, dimension (N) - On entry, if JPVT(i) .ne. 0, the i-th column of A is an - initial column, otherwise it is a free column. Before - the QR factorization of A, all initial columns are - permuted to the leading positions; only the remaining - free columns are moved as a result of column pivoting - during the factorization. - On exit, if JPVT(i) = k, then the i-th column of A*P - was the k-th column of A. - - RCOND (input) DOUBLE PRECISION - RCOND is used to determine the effective rank of A, which - is defined as the order of the largest leading triangular - submatrix R11 in the QR factorization with pivoting of A, - whose estimated condition number < 1/RCOND. - - RANK (output) INTEGER - The effective rank of A, i.e., the order of the submatrix - R11. This is the same as the order of the submatrix T11 - in the complete orthogonal factorization of A. - - WORK (workspace) DOUBLE PRECISION array, dimension - (max( min(M,N)+3*N, 2*min(M,N)+NRHS )), - - INFO (output) INTEGER - = 0: successful exit - < 0: if INFO = -i, the i-th argument had an illegal value - - ===================================================================== - - - Parameter adjustments */ - /* Table of constant values */ - static integer c__0 = 0; - static doublereal c_b13 = 0.; - static integer c__2 = 2; - static integer c__1 = 1; - static doublereal c_b36 = 1.; - - /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; - doublereal d__1; - /* Local variables */ - static doublereal anrm, bnrm, smin, smax; - static integer i__, j, k, iascl, ibscl, ismin, ismax; - static doublereal c1, c2; - extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *, - integer *, integer *, doublereal *, doublereal *, integer *, - doublereal *, integer *), dlaic1_( - integer *, integer *, doublereal *, doublereal *, doublereal *, - doublereal *, doublereal *, doublereal *, doublereal *); - static doublereal s1, s2, t1, t2; - extern /* Subroutine */ int dorm2r_(char *, char *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - integer *, doublereal *, integer *), dlabad_( - doublereal *, doublereal *); - extern doublereal dlamch_(char *), dlange_(char *, integer *, - integer *, doublereal *, integer *, doublereal *); - static integer mn; - extern /* Subroutine */ int dlascl_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, integer *, doublereal *, - integer *, integer *), dgeqpf_(integer *, integer *, - doublereal *, integer *, integer *, doublereal *, doublereal *, - integer *), dlaset_(char *, integer *, integer *, doublereal *, - doublereal *, doublereal *, integer *), xerbla_(char *, - integer *); - static doublereal bignum; - extern /* Subroutine */ int dlatzm_(char *, integer *, integer *, - doublereal *, integer *, doublereal *, doublereal *, doublereal *, - integer *, doublereal *); - static doublereal sminpr, smaxpr, smlnum; - extern /* Subroutine */ int dtzrqf_(integer *, integer *, doublereal *, - integer *, doublereal *, integer *); -#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] -#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] - - - a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; - a -= a_offset; - b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; - b -= b_offset; - --jpvt; - --work; - - /* Function Body */ - mn = min(*m,*n); - ismin = mn + 1; - ismax = (mn << 1) + 1; - -/* Test the input arguments. */ - - *info = 0; - if (*m < 0) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*nrhs < 0) { - *info = -3; - } else if (*lda < max(1,*m)) { - *info = -5; - } else /* if(complicated condition) */ { -/* Computing MAX */ - i__1 = max(1,*m); - if (*ldb < max(i__1,*n)) { - *info = -7; - } - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("DGELSX", &i__1); - return 0; - } - -/* Quick return if possible - - Computing MIN */ - i__1 = min(*m,*n); - if (min(i__1,*nrhs) == 0) { - *rank = 0; - return 0; - } - -/* Get machine parameters */ - - smlnum = dlamch_("S") / dlamch_("P"); - bignum = 1. / smlnum; - dlabad_(&smlnum, &bignum); - -/* Scale A, B if max elements outside range [SMLNUM,BIGNUM] */ - - anrm = dlange_("M", m, n, &a[a_offset], lda, &work[1]); - iascl = 0; - if (anrm > 0. && anrm < smlnum) { - -/* Scale matrix norm up to SMLNUM */ - - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, - info); - iascl = 1; - } else if (anrm > bignum) { - -/* Scale matrix norm down to BIGNUM */ - - dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, - info); - iascl = 2; - } else if (anrm == 0.) { - -/* Matrix all zero. Return zero solution. */ - - i__1 = max(*m,*n); - dlaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb); - *rank = 0; - goto L100; - } - - bnrm = dlange_("M", m, nrhs, &b[b_offset], ldb, &work[1]); - ibscl = 0; - if (bnrm > 0. && bnrm < smlnum) { - -/* Scale matrix norm up to SMLNUM */ - - dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 1; - } else if (bnrm > bignum) { - -/* Scale matrix norm down to BIGNUM */ - - dlascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 2; - } - -/* Compute QR factorization with column pivoting of A: - A * P = Q * R */ - - dgeqpf_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], info); - -/* workspace 3*N. Details of Householder rotations stored - in WORK(1:MN). - - Determine RANK using incremental condition estimation */ - - work[ismin] = 1.; - work[ismax] = 1.; - smax = (d__1 = a_ref(1, 1), abs(d__1)); - smin = smax; - if ((d__1 = a_ref(1, 1), abs(d__1)) == 0.) { - *rank = 0; - i__1 = max(*m,*n); - dlaset_("F", &i__1, nrhs, &c_b13, &c_b13, &b[b_offset], ldb); - goto L100; - } else { - *rank = 1; - } - -L10: - if (*rank < mn) { - i__ = *rank + 1; - dlaic1_(&c__2, rank, &work[ismin], &smin, &a_ref(1, i__), &a_ref(i__, - i__), &sminpr, &s1, &c1); - dlaic1_(&c__1, rank, &work[ismax], &smax, &a_ref(1, i__), &a_ref(i__, - i__), &smaxpr, &s2, &c2); - - if (smaxpr * *rcond <= sminpr) { - i__1 = *rank; - for (i__ = 1; i__ <= i__1; ++i__) { - work[ismin + i__ - 1] = s1 * work[ismin + i__ - 1]; - work[ismax + i__ - 1] = s2 * work[ismax + i__ - 1]; -/* L20: */ - } - work[ismin + *rank] = c1; - work[ismax + *rank] = c2; - smin = sminpr; - smax = smaxpr; - ++(*rank); - goto L10; - } - } - -/* Logically partition R = [ R11 R12 ] - [ 0 R22 ] - where R11 = R(1:RANK,1:RANK) - - [R11,R12] = [ T11, 0 ] * Y */ - - if (*rank < *n) { - dtzrqf_(rank, n, &a[a_offset], lda, &work[mn + 1], info); - } - -/* Details of Householder rotations stored in WORK(MN+1:2*MN) - - B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */ - - dorm2r_("Left", "Transpose", m, nrhs, &mn, &a[a_offset], lda, &work[1], & - b[b_offset], ldb, &work[(mn << 1) + 1], info); - -/* workspace NRHS - - B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */ - - dtrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b36, & - a[a_offset], lda, &b[b_offset], ldb); - - i__1 = *n; - for (i__ = *rank + 1; i__ <= i__1; ++i__) { - i__2 = *nrhs; - for (j = 1; j <= i__2; ++j) { - b_ref(i__, j) = 0.; -/* L30: */ - } -/* L40: */ - } - -/* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */ - - if (*rank < *n) { - i__1 = *rank; - for (i__ = 1; i__ <= i__1; ++i__) { - i__2 = *n - *rank + 1; - dlatzm_("Left", &i__2, nrhs, &a_ref(i__, *rank + 1), lda, &work[ - mn + i__], &b_ref(i__, 1), &b_ref(*rank + 1, 1), ldb, & - work[(mn << 1) + 1]); -/* L50: */ - } - } - -/* workspace NRHS - - B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */ - - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - work[(mn << 1) + i__] = 1.; -/* L60: */ - } - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - if (work[(mn << 1) + i__] == 1.) { - if (jpvt[i__] != i__) { - k = i__; - t1 = b_ref(k, j); - t2 = b_ref(jpvt[k], j); -L70: - b_ref(jpvt[k], j) = t1; - work[(mn << 1) + k] = 0.; - t1 = t2; - k = jpvt[k]; - t2 = b_ref(jpvt[k], j); - if (jpvt[k] != i__) { - goto L70; - } - b_ref(i__, j) = t1; - work[(mn << 1) + k] = 0.; - } - } -/* L80: */ - } -/* L90: */ - } - -/* Undo scaling */ - - if (iascl == 1) { - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset], - lda, info); - } else if (iascl == 2) { - dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset], - lda, info); - } - if (ibscl == 1) { - dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } else if (ibscl == 2) { - dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } - -L100: - - return 0; - -/* End of DGELSX */ - -} /* dgelsx_ */ - -#undef b_ref -#undef a_ref - - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_lapack/dgelsy.c b/ext/f2c_lapack/dgelsy.c deleted file mode 100644 index 13be12be7..000000000 --- a/ext/f2c_lapack/dgelsy.c +++ /dev/null @@ -1,475 +0,0 @@ -#include "blaswrap.h" -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -/* Subroutine */ int dgelsy_(integer *m, integer *n, integer *nrhs, - doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * - jpvt, doublereal *rcond, integer *rank, doublereal *work, integer * - lwork, integer *info) -{ -/* -- LAPACK driver routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - June 30, 1999 - - - Purpose - ======= - - DGELSY computes the minimum-norm solution to a real linear least - squares problem: - minimize || A * X - B || - using a complete orthogonal factorization of A. A is an M-by-N - matrix which may be rank-deficient. - - Several right hand side vectors b and solution vectors x can be - handled in a single call; they are stored as the columns of the - M-by-NRHS right hand side matrix B and the N-by-NRHS solution - matrix X. - - The routine first computes a QR factorization with column pivoting: - A * P = Q * [ R11 R12 ] - [ 0 R22 ] - with R11 defined as the largest leading submatrix whose estimated - condition number is less than 1/RCOND. The order of R11, RANK, - is the effective rank of A. - - Then, R22 is considered to be negligible, and R12 is annihilated - by orthogonal transformations from the right, arriving at the - complete orthogonal factorization: - A * P = Q * [ T11 0 ] * Z - [ 0 0 ] - The minimum-norm solution is then - X = P * Z' [ inv(T11)*Q1'*B ] - [ 0 ] - where Q1 consists of the first RANK columns of Q. - - This routine is basically identical to the original xGELSX except - three differences: - o The call to the subroutine xGEQPF has been substituted by the - the call to the subroutine xGEQP3. This subroutine is a Blas-3 - version of the QR factorization with column pivoting. - o Matrix B (the right hand side) is updated with Blas-3. - o The permutation of matrix B (the right hand side) is faster and - more simple. - - Arguments - ========= - - M (input) INTEGER - The number of rows of the matrix A. M >= 0. - - N (input) INTEGER - The number of columns of the matrix A. N >= 0. - - NRHS (input) INTEGER - The number of right hand sides, i.e., the number of - columns of matrices B and X. NRHS >= 0. - - A (input/output) DOUBLE PRECISION array, dimension (LDA,N) - On entry, the M-by-N matrix A. - On exit, A has been overwritten by details of its - complete orthogonal factorization. - - LDA (input) INTEGER - The leading dimension of the array A. LDA >= max(1,M). - - B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS) - On entry, the M-by-NRHS right hand side matrix B. - On exit, the N-by-NRHS solution matrix X. - - LDB (input) INTEGER - The leading dimension of the array B. LDB >= max(1,M,N). - - JPVT (input/output) INTEGER array, dimension (N) - On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted - to the front of AP, otherwise column i is a free column. - On exit, if JPVT(i) = k, then the i-th column of AP - was the k-th column of A. - - RCOND (input) DOUBLE PRECISION - RCOND is used to determine the effective rank of A, which - is defined as the order of the largest leading triangular - submatrix R11 in the QR factorization with pivoting of A, - whose estimated condition number < 1/RCOND. - - RANK (output) INTEGER - The effective rank of A, i.e., the order of the submatrix - R11. This is the same as the order of the submatrix T11 - in the complete orthogonal factorization of A. - - WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) - On exit, if INFO = 0, WORK(1) returns the optimal LWORK. - - LWORK (input) INTEGER - The dimension of the array WORK. - The unblocked strategy requires that: - LWORK >= MAX( MN+3*N+1, 2*MN+NRHS ), - where MN = min( M, N ). - The block algorithm requires that: - LWORK >= MAX( MN+2*N+NB*(N+1), 2*MN+NB*NRHS ), - where NB is an upper bound on the blocksize returned - by ILAENV for the routines DGEQP3, DTZRZF, STZRQF, DORMQR, - and DORMRZ. - - If LWORK = -1, then a workspace query is assumed; the routine - only calculates the optimal size of the WORK array, returns - this value as the first entry of the WORK array, and no error - message related to LWORK is issued by XERBLA. - - INFO (output) INTEGER - = 0: successful exit - < 0: If INFO = -i, the i-th argument had an illegal value. - - Further Details - =============== - - Based on contributions by - A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA - E. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain - G. Quintana-Orti, Depto. de Informatica, Universidad Jaime I, Spain - - ===================================================================== - - - Parameter adjustments */ - /* Table of constant values */ - static integer c__1 = 1; - static integer c_n1 = -1; - static integer c__0 = 0; - static doublereal c_b31 = 0.; - static integer c__2 = 2; - static doublereal c_b54 = 1.; - - /* System generated locals */ - integer a_dim1, a_offset, b_dim1, b_offset, i__1, i__2; - doublereal d__1, d__2; - /* Local variables */ - static doublereal anrm, bnrm, smin, smax; - static integer i__, j, iascl, ibscl; - extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, - doublereal *, integer *); - static integer ismin, ismax; - static doublereal c1, c2; - extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *, - integer *, integer *, doublereal *, doublereal *, integer *, - doublereal *, integer *), dlaic1_( - integer *, integer *, doublereal *, doublereal *, doublereal *, - doublereal *, doublereal *, doublereal *, doublereal *); - static doublereal wsize, s1, s2; - extern /* Subroutine */ int dgeqp3_(integer *, integer *, doublereal *, - integer *, integer *, doublereal *, doublereal *, integer *, - integer *), dlabad_(doublereal *, doublereal *); - static integer nb; - extern doublereal dlamch_(char *), dlange_(char *, integer *, - integer *, doublereal *, integer *, doublereal *); - static integer mn; - extern /* Subroutine */ int dlascl_(char *, integer *, integer *, - doublereal *, doublereal *, integer *, integer *, doublereal *, - integer *, integer *), dlaset_(char *, integer *, integer - *, doublereal *, doublereal *, doublereal *, integer *), - xerbla_(char *, integer *); - extern integer ilaenv_(integer *, char *, char *, integer *, integer *, - integer *, integer *, ftnlen, ftnlen); - static doublereal bignum; - static integer nb1, nb2, nb3, nb4; - extern /* Subroutine */ int dormqr_(char *, char *, integer *, integer *, - integer *, doublereal *, integer *, doublereal *, doublereal *, - integer *, doublereal *, integer *, integer *); - static doublereal sminpr, smaxpr, smlnum; - extern /* Subroutine */ int dormrz_(char *, char *, integer *, integer *, - integer *, integer *, doublereal *, integer *, doublereal *, - doublereal *, integer *, doublereal *, integer *, integer *); - static integer lwkopt; - static logical lquery; - extern /* Subroutine */ int dtzrzf_(integer *, integer *, doublereal *, - integer *, doublereal *, doublereal *, integer *, integer *); -#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] -#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] - - - a_dim1 = *lda; - a_offset = 1 + a_dim1 * 1; - a -= a_offset; - b_dim1 = *ldb; - b_offset = 1 + b_dim1 * 1; - b -= b_offset; - --jpvt; - --work; - - /* Function Body */ - mn = min(*m,*n); - ismin = mn + 1; - ismax = (mn << 1) + 1; - -/* Test the input arguments. */ - - *info = 0; - nb1 = ilaenv_(&c__1, "DGEQRF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, ( - ftnlen)1); - nb2 = ilaenv_(&c__1, "DGERQF", " ", m, n, &c_n1, &c_n1, (ftnlen)6, ( - ftnlen)1); - nb3 = ilaenv_(&c__1, "DORMQR", " ", m, n, nrhs, &c_n1, (ftnlen)6, (ftnlen) - 1); - nb4 = ilaenv_(&c__1, "DORMRQ", " ", m, n, nrhs, &c_n1, (ftnlen)6, (ftnlen) - 1); -/* Computing MAX */ - i__1 = max(nb1,nb2), i__1 = max(i__1,nb3); - nb = max(i__1,nb4); -/* Computing MAX */ - i__1 = 1, i__2 = mn + (*n << 1) + nb * (*n + 1), i__1 = max(i__1,i__2), - i__2 = (mn << 1) + nb * *nrhs; - lwkopt = max(i__1,i__2); - work[1] = (doublereal) lwkopt; - lquery = *lwork == -1; - if (*m < 0) { - *info = -1; - } else if (*n < 0) { - *info = -2; - } else if (*nrhs < 0) { - *info = -3; - } else if (*lda < max(1,*m)) { - *info = -5; - } else /* if(complicated condition) */ { -/* Computing MAX */ - i__1 = max(1,*m); - if (*ldb < max(i__1,*n)) { - *info = -7; - } else /* if(complicated condition) */ { -/* Computing MAX */ - i__1 = 1, i__2 = mn + *n * 3 + 1, i__1 = max(i__1,i__2), i__2 = ( - mn << 1) + *nrhs; - if (*lwork < max(i__1,i__2) && ! lquery) { - *info = -12; - } - } - } - - if (*info != 0) { - i__1 = -(*info); - xerbla_("DGELSY", &i__1); - return 0; - } else if (lquery) { - return 0; - } - -/* Quick return if possible - - Computing MIN */ - i__1 = min(*m,*n); - if (min(i__1,*nrhs) == 0) { - *rank = 0; - return 0; - } - -/* Get machine parameters */ - - smlnum = dlamch_("S") / dlamch_("P"); - bignum = 1. / smlnum; - dlabad_(&smlnum, &bignum); - -/* Scale A, B if max entries outside range [SMLNUM,BIGNUM] */ - - anrm = dlange_("M", m, n, &a[a_offset], lda, &work[1]); - iascl = 0; - if (anrm > 0. && anrm < smlnum) { - -/* Scale matrix norm up to SMLNUM */ - - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, m, n, &a[a_offset], lda, - info); - iascl = 1; - } else if (anrm > bignum) { - -/* Scale matrix norm down to BIGNUM */ - - dlascl_("G", &c__0, &c__0, &anrm, &bignum, m, n, &a[a_offset], lda, - info); - iascl = 2; - } else if (anrm == 0.) { - -/* Matrix all zero. Return zero solution. */ - - i__1 = max(*m,*n); - dlaset_("F", &i__1, nrhs, &c_b31, &c_b31, &b[b_offset], ldb); - *rank = 0; - goto L70; - } - - bnrm = dlange_("M", m, nrhs, &b[b_offset], ldb, &work[1]); - ibscl = 0; - if (bnrm > 0. && bnrm < smlnum) { - -/* Scale matrix norm up to SMLNUM */ - - dlascl_("G", &c__0, &c__0, &bnrm, &smlnum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 1; - } else if (bnrm > bignum) { - -/* Scale matrix norm down to BIGNUM */ - - dlascl_("G", &c__0, &c__0, &bnrm, &bignum, m, nrhs, &b[b_offset], ldb, - info); - ibscl = 2; - } - -/* Compute QR factorization with column pivoting of A: - A * P = Q * R */ - - i__1 = *lwork - mn; - dgeqp3_(m, n, &a[a_offset], lda, &jpvt[1], &work[1], &work[mn + 1], &i__1, - info); - wsize = mn + work[mn + 1]; - -/* workspace: MN+2*N+NB*(N+1). - Details of Householder rotations stored in WORK(1:MN). - - Determine RANK using incremental condition estimation */ - - work[ismin] = 1.; - work[ismax] = 1.; - smax = (d__1 = a_ref(1, 1), abs(d__1)); - smin = smax; - if ((d__1 = a_ref(1, 1), abs(d__1)) == 0.) { - *rank = 0; - i__1 = max(*m,*n); - dlaset_("F", &i__1, nrhs, &c_b31, &c_b31, &b[b_offset], ldb); - goto L70; - } else { - *rank = 1; - } - -L10: - if (*rank < mn) { - i__ = *rank + 1; - dlaic1_(&c__2, rank, &work[ismin], &smin, &a_ref(1, i__), &a_ref(i__, - i__), &sminpr, &s1, &c1); - dlaic1_(&c__1, rank, &work[ismax], &smax, &a_ref(1, i__), &a_ref(i__, - i__), &smaxpr, &s2, &c2); - - if (smaxpr * *rcond <= sminpr) { - i__1 = *rank; - for (i__ = 1; i__ <= i__1; ++i__) { - work[ismin + i__ - 1] = s1 * work[ismin + i__ - 1]; - work[ismax + i__ - 1] = s2 * work[ismax + i__ - 1]; -/* L20: */ - } - work[ismin + *rank] = c1; - work[ismax + *rank] = c2; - smin = sminpr; - smax = smaxpr; - ++(*rank); - goto L10; - } - } - -/* workspace: 3*MN. - - Logically partition R = [ R11 R12 ] - [ 0 R22 ] - where R11 = R(1:RANK,1:RANK) - - [R11,R12] = [ T11, 0 ] * Y */ - - if (*rank < *n) { - i__1 = *lwork - (mn << 1); - dtzrzf_(rank, n, &a[a_offset], lda, &work[mn + 1], &work[(mn << 1) + - 1], &i__1, info); - } - -/* workspace: 2*MN. - Details of Householder rotations stored in WORK(MN+1:2*MN) - - B(1:M,1:NRHS) := Q' * B(1:M,1:NRHS) */ - - i__1 = *lwork - (mn << 1); - dormqr_("Left", "Transpose", m, nrhs, &mn, &a[a_offset], lda, &work[1], & - b[b_offset], ldb, &work[(mn << 1) + 1], &i__1, info); -/* Computing MAX */ - d__1 = wsize, d__2 = (mn << 1) + work[(mn << 1) + 1]; - wsize = max(d__1,d__2); - -/* workspace: 2*MN+NB*NRHS. - - B(1:RANK,1:NRHS) := inv(T11) * B(1:RANK,1:NRHS) */ - - dtrsm_("Left", "Upper", "No transpose", "Non-unit", rank, nrhs, &c_b54, & - a[a_offset], lda, &b[b_offset], ldb); - - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = *rank + 1; i__ <= i__2; ++i__) { - b_ref(i__, j) = 0.; -/* L30: */ - } -/* L40: */ - } - -/* B(1:N,1:NRHS) := Y' * B(1:N,1:NRHS) */ - - if (*rank < *n) { - i__1 = *n - *rank; - i__2 = *lwork - (mn << 1); - dormrz_("Left", "Transpose", n, nrhs, rank, &i__1, &a[a_offset], lda, - &work[mn + 1], &b[b_offset], ldb, &work[(mn << 1) + 1], &i__2, - info); - } - -/* workspace: 2*MN+NRHS. - - B(1:N,1:NRHS) := P * B(1:N,1:NRHS) */ - - i__1 = *nrhs; - for (j = 1; j <= i__1; ++j) { - i__2 = *n; - for (i__ = 1; i__ <= i__2; ++i__) { - work[jpvt[i__]] = b_ref(i__, j); -/* L50: */ - } - dcopy_(n, &work[1], &c__1, &b_ref(1, j), &c__1); -/* L60: */ - } - -/* workspace: N. - - Undo scaling */ - - if (iascl == 1) { - dlascl_("G", &c__0, &c__0, &anrm, &smlnum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("U", &c__0, &c__0, &smlnum, &anrm, rank, rank, &a[a_offset], - lda, info); - } else if (iascl == 2) { - dlascl_("G", &c__0, &c__0, &anrm, &bignum, n, nrhs, &b[b_offset], ldb, - info); - dlascl_("U", &c__0, &c__0, &bignum, &anrm, rank, rank, &a[a_offset], - lda, info); - } - if (ibscl == 1) { - dlascl_("G", &c__0, &c__0, &smlnum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } else if (ibscl == 2) { - dlascl_("G", &c__0, &c__0, &bignum, &bnrm, n, nrhs, &b[b_offset], ldb, - info); - } - -L70: - work[1] = (doublereal) lwkopt; - - return 0; - -/* End of DGELSY */ - -} /* dgelsy_ */ - -#undef b_ref -#undef a_ref - - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_lapack/lsame.c b/ext/f2c_lapack/lsame.c deleted file mode 100644 index ba8740b44..000000000 --- a/ext/f2c_lapack/lsame.c +++ /dev/null @@ -1,107 +0,0 @@ -#ifdef _cpluscplus -extern "C" { -#endif -#include "f2c.h" - -logical lsame_(char *ca, char *cb) -{ -/* -- LAPACK auxiliary routine (version 3.0) -- - Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., - Courant Institute, Argonne National Lab, and Rice University - September 30, 1994 - - - Purpose - ======= - - LSAME returns .TRUE. if CA is the same letter as CB regardless of - case. - - Arguments - ========= - - CA (input) CHARACTER*1 - CB (input) CHARACTER*1 - CA and CB specify the single characters to be compared. - - ===================================================================== - - - - Test if the characters are equal */ - /* System generated locals */ - logical ret_val; - /* Local variables */ - static integer inta, intb, zcode; - - - ret_val = *(unsigned char *)ca == *(unsigned char *)cb; - if (ret_val) { - return ret_val; - } - -/* Now test for equivalence if both characters are alphabetic. */ - - zcode = 'Z'; - -/* Use 'Z' rather than 'A' so that ASCII can be detected on Prime - machines, on which ICHAR returns a value with bit 8 set. - ICHAR('A') on Prime machines returns 193 which is the same as - ICHAR('A') on an EBCDIC machine. */ - - inta = *(unsigned char *)ca; - intb = *(unsigned char *)cb; - - if (zcode == 90 || zcode == 122) { - -/* ASCII is assumed - ZCODE is the ASCII code of either lower o -r - upper case 'Z'. */ - - if (inta >= 97 && inta <= 122) { - inta += -32; - } - if (intb >= 97 && intb <= 122) { - intb += -32; - } - - } else if (zcode == 233 || zcode == 169) { - -/* EBCDIC is assumed - ZCODE is the EBCDIC code of either lower - or - upper case 'Z'. */ - - if (inta >= 129 && inta <= 137 || inta >= 145 && inta <= 153 || inta - >= 162 && inta <= 169) { - inta += 64; - } - if (intb >= 129 && intb <= 137 || intb >= 145 && intb <= 153 || intb - >= 162 && intb <= 169) { - intb += 64; - } - - } else if (zcode == 218 || zcode == 250) { - -/* ASCII is assumed, on Prime machines - ZCODE is the ASCII cod -e - plus 128 of either lower or upper case 'Z'. */ - - if (inta >= 225 && inta <= 250) { - inta += -32; - } - if (intb >= 225 && intb <= 250) { - intb += -32; - } - } - ret_val = inta == intb; - -/* RETURN - - End of LSAME */ - - return ret_val; -} /* lsame_ */ - -#ifdef _cpluscplus -} -#endif diff --git a/ext/f2c_libs/derf_.c b/ext/f2c_libs/derf_.c deleted file mode 100644 index d935d3152..000000000 --- a/ext/f2c_libs/derf_.c +++ /dev/null @@ -1,18 +0,0 @@ -#include "f2c.h" -#ifdef __cplusplus -extern "C" { -#endif - -#ifdef KR_headers -double erf(); -double derf_(x) doublereal *x; -#else -extern double erf(double); -double derf_(doublereal *x) -#endif -{ -return( erf(*x) ); -} -#ifdef __cplusplus -} -#endif diff --git a/ext/f2c_libs/derfc_.c b/ext/f2c_libs/derfc_.c deleted file mode 100644 index 18f5c619b..000000000 --- a/ext/f2c_libs/derfc_.c +++ /dev/null @@ -1,20 +0,0 @@ -#include "f2c.h" -#ifdef __cplusplus -extern "C" { -#endif - -#ifdef KR_headers -extern double erfc(); - -double derfc_(x) doublereal *x; -#else -extern double erfc(double); - -double derfc_(doublereal *x) -#endif -{ -return( erfc(*x) ); -} -#ifdef __cplusplus -} -#endif diff --git a/ext/f2c_libs/erf_.c b/ext/f2c_libs/erf_.c deleted file mode 100644 index 532fec61c..000000000 --- a/ext/f2c_libs/erf_.c +++ /dev/null @@ -1,22 +0,0 @@ -#include "f2c.h" -#ifdef __cplusplus -extern "C" { -#endif - -#ifndef REAL -#define REAL double -#endif - -#ifdef KR_headers -double erf(); -REAL erf_(x) real *x; -#else -extern double erf(double); -REAL erf_(real *x) -#endif -{ -return( erf((double)*x) ); -} -#ifdef __cplusplus -} -#endif diff --git a/ext/f2c_libs/erfc_.c b/ext/f2c_libs/erfc_.c deleted file mode 100644 index 6f6c9f106..000000000 --- a/ext/f2c_libs/erfc_.c +++ /dev/null @@ -1,22 +0,0 @@ -#include "f2c.h" -#ifdef __cplusplus -extern "C" { -#endif - -#ifndef REAL -#define REAL double -#endif - -#ifdef KR_headers -double erfc(); -REAL erfc_(x) real *x; -#else -extern double erfc(double); -REAL erfc_(real *x) -#endif -{ -return( erfc((double)*x) ); -} -#ifdef __cplusplus -} -#endif diff --git a/ext/f2c_libs/getarg_.c b/ext/f2c_libs/getarg_.c deleted file mode 100644 index 9006ad201..000000000 --- a/ext/f2c_libs/getarg_.c +++ /dev/null @@ -1,34 +0,0 @@ -#include "f2c.h" -#ifdef __cplusplus -extern "C" { -#endif - -/* - * subroutine getarg(k, c) - * returns the kth unix command argument in fortran character - * variable argument c -*/ - -#ifdef KR_headers -VOID getarg_(n, s, ls) ftnint *n; register char *s; ftnlen ls; -#else -void getarg_(ftnint *n, register char *s, ftnlen ls) -#endif -{ -extern int xargc; -extern char **xargv; -register char *t; -register int i; - -if(*n>=0 && *n