sgesc2.c
Go to the documentation of this file.
00001 /* sgesc2.f -- translated by f2c (version 20061008).
00002    You must link the resulting object file with libf2c:
00003         on Microsoft Windows system, link with libf2c.lib;
00004         on Linux or Unix systems, link with .../path/to/libf2c.a -lm
00005         or, if you install libf2c.a in a standard place, with -lf2c -lm
00006         -- in that order, at the end of the command line, as in
00007                 cc *.o -lf2c -lm
00008         Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
00009 
00010                 http://www.netlib.org/f2c/libf2c.zip
00011 */
00012 
00013 #include "f2c.h"
00014 #include "blaswrap.h"
00015 
00016 /* Table of constant values */
00017 
00018 static integer c__1 = 1;
00019 static integer c_n1 = -1;
00020 
00021 /* Subroutine */ int sgesc2_(integer *n, real *a, integer *lda, real *rhs, 
00022         integer *ipiv, integer *jpiv, real *scale)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, i__1, i__2;
00026     real r__1, r__2;
00027 
00028     /* Local variables */
00029     integer i__, j;
00030     real eps, temp;
00031     extern /* Subroutine */ int sscal_(integer *, real *, real *, integer *), 
00032             slabad_(real *, real *);
00033     extern doublereal slamch_(char *);
00034     real bignum;
00035     extern integer isamax_(integer *, real *, integer *);
00036     extern /* Subroutine */ int slaswp_(integer *, real *, integer *, integer 
00037             *, integer *, integer *, integer *);
00038     real smlnum;
00039 
00040 
00041 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00042 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00043 /*     November 2006 */
00044 
00045 /*     .. Scalar Arguments .. */
00046 /*     .. */
00047 /*     .. Array Arguments .. */
00048 /*     .. */
00049 
00050 /*  Purpose */
00051 /*  ======= */
00052 
00053 /*  SGESC2 solves a system of linear equations */
00054 
00055 /*            A * X = scale* RHS */
00056 
00057 /*  with a general N-by-N matrix A using the LU factorization with */
00058 /*  complete pivoting computed by SGETC2. */
00059 
00060 /*  Arguments */
00061 /*  ========= */
00062 
00063 /*  N       (input) INTEGER */
00064 /*          The order of the matrix A. */
00065 
00066 /*  A       (input) REAL array, dimension (LDA,N) */
00067 /*          On entry, the  LU part of the factorization of the n-by-n */
00068 /*          matrix A computed by SGETC2:  A = P * L * U * Q */
00069 
00070 /*  LDA     (input) INTEGER */
00071 /*          The leading dimension of the array A.  LDA >= max(1, N). */
00072 
00073 /*  RHS     (input/output) REAL array, dimension (N). */
00074 /*          On entry, the right hand side vector b. */
00075 /*          On exit, the solution vector X. */
00076 
00077 /*  IPIV    (input) INTEGER array, dimension (N). */
00078 /*          The pivot indices; for 1 <= i <= N, row i of the */
00079 /*          matrix has been interchanged with row IPIV(i). */
00080 
00081 /*  JPIV    (input) INTEGER array, dimension (N). */
00082 /*          The pivot indices; for 1 <= j <= N, column j of the */
00083 /*          matrix has been interchanged with column JPIV(j). */
00084 
00085 /*  SCALE    (output) REAL */
00086 /*           On exit, SCALE contains the scale factor. SCALE is chosen */
00087 /*           0 <= SCALE <= 1 to prevent owerflow in the solution. */
00088 
00089 /*  Further Details */
00090 /*  =============== */
00091 
00092 /*  Based on contributions by */
00093 /*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
00094 /*     Umea University, S-901 87 Umea, Sweden. */
00095 
00096 /*  ===================================================================== */
00097 
00098 /*     .. Parameters .. */
00099 /*     .. */
00100 /*     .. Local Scalars .. */
00101 /*     .. */
00102 /*     .. External Subroutines .. */
00103 /*     .. */
00104 /*     .. External Functions .. */
00105 /*     .. */
00106 /*     .. Intrinsic Functions .. */
00107 /*     .. */
00108 /*     .. Executable Statements .. */
00109 
00110 /*      Set constant to control owerflow */
00111 
00112     /* Parameter adjustments */
00113     a_dim1 = *lda;
00114     a_offset = 1 + a_dim1;
00115     a -= a_offset;
00116     --rhs;
00117     --ipiv;
00118     --jpiv;
00119 
00120     /* Function Body */
00121     eps = slamch_("P");
00122     smlnum = slamch_("S") / eps;
00123     bignum = 1.f / smlnum;
00124     slabad_(&smlnum, &bignum);
00125 
00126 /*     Apply permutations IPIV to RHS */
00127 
00128     i__1 = *n - 1;
00129     slaswp_(&c__1, &rhs[1], lda, &c__1, &i__1, &ipiv[1], &c__1);
00130 
00131 /*     Solve for L part */
00132 
00133     i__1 = *n - 1;
00134     for (i__ = 1; i__ <= i__1; ++i__) {
00135         i__2 = *n;
00136         for (j = i__ + 1; j <= i__2; ++j) {
00137             rhs[j] -= a[j + i__ * a_dim1] * rhs[i__];
00138 /* L10: */
00139         }
00140 /* L20: */
00141     }
00142 
00143 /*     Solve for U part */
00144 
00145     *scale = 1.f;
00146 
00147 /*     Check for scaling */
00148 
00149     i__ = isamax_(n, &rhs[1], &c__1);
00150     if (smlnum * 2.f * (r__1 = rhs[i__], dabs(r__1)) > (r__2 = a[*n + *n * 
00151             a_dim1], dabs(r__2))) {
00152         temp = .5f / (r__1 = rhs[i__], dabs(r__1));
00153         sscal_(n, &temp, &rhs[1], &c__1);
00154         *scale *= temp;
00155     }
00156 
00157     for (i__ = *n; i__ >= 1; --i__) {
00158         temp = 1.f / a[i__ + i__ * a_dim1];
00159         rhs[i__] *= temp;
00160         i__1 = *n;
00161         for (j = i__ + 1; j <= i__1; ++j) {
00162             rhs[i__] -= rhs[j] * (a[i__ + j * a_dim1] * temp);
00163 /* L30: */
00164         }
00165 /* L40: */
00166     }
00167 
00168 /*     Apply permutations JPIV to the solution (RHS) */
00169 
00170     i__1 = *n - 1;
00171     slaswp_(&c__1, &rhs[1], lda, &c__1, &i__1, &jpiv[1], &c_n1);
00172     return 0;
00173 
00174 /*     End of SGESC2 */
00175 
00176 } /* sgesc2_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:56:07