dgetc2.c
Go to the documentation of this file.
00001 /* dgetc2.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 doublereal c_b10 = -1.;
00020 
00021 /* Subroutine */ int dgetc2_(integer *n, doublereal *a, integer *lda, integer 
00022         *ipiv, integer *jpiv, integer *info)
00023 {
00024     /* System generated locals */
00025     integer a_dim1, a_offset, i__1, i__2, i__3;
00026     doublereal d__1;
00027 
00028     /* Local variables */
00029     integer i__, j, ip, jp;
00030     doublereal eps;
00031     integer ipv, jpv;
00032     extern /* Subroutine */ int dger_(integer *, integer *, doublereal *, 
00033             doublereal *, integer *, doublereal *, integer *, doublereal *, 
00034             integer *);
00035     doublereal smin, xmax;
00036     extern /* Subroutine */ int dswap_(integer *, doublereal *, integer *, 
00037             doublereal *, integer *), dlabad_(doublereal *, doublereal *);
00038     extern doublereal dlamch_(char *);
00039     doublereal bignum, smlnum;
00040 
00041 
00042 /*  -- LAPACK auxiliary routine (version 3.2) -- */
00043 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
00044 /*     November 2006 */
00045 
00046 /*     .. Scalar Arguments .. */
00047 /*     .. */
00048 /*     .. Array Arguments .. */
00049 /*     .. */
00050 
00051 /*  Purpose */
00052 /*  ======= */
00053 
00054 /*  DGETC2 computes an LU factorization with complete pivoting of the */
00055 /*  n-by-n matrix A. The factorization has the form A = P * L * U * Q, */
00056 /*  where P and Q are permutation matrices, L is lower triangular with */
00057 /*  unit diagonal elements and U is upper triangular. */
00058 
00059 /*  This is the Level 2 BLAS algorithm. */
00060 
00061 /*  Arguments */
00062 /*  ========= */
00063 
00064 /*  N       (input) INTEGER */
00065 /*          The order of the matrix A. N >= 0. */
00066 
00067 /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N) */
00068 /*          On entry, the n-by-n matrix A to be factored. */
00069 /*          On exit, the factors L and U from the factorization */
00070 /*          A = P*L*U*Q; the unit diagonal elements of L are not stored. */
00071 /*          If U(k, k) appears to be less than SMIN, U(k, k) is given the */
00072 /*          value of SMIN, i.e., giving a nonsingular perturbed system. */
00073 
00074 /*  LDA     (input) INTEGER */
00075 /*          The leading dimension of the array A.  LDA >= max(1,N). */
00076 
00077 /*  IPIV    (output) 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    (output) 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 /*  INFO    (output) INTEGER */
00086 /*           = 0: successful exit */
00087 /*           > 0: if INFO = k, U(k, k) is likely to produce owerflow if */
00088 /*                we try to solve for x in Ax = b. So U is perturbed to */
00089 /*                avoid the overflow. */
00090 
00091 /*  Further Details */
00092 /*  =============== */
00093 
00094 /*  Based on contributions by */
00095 /*     Bo Kagstrom and Peter Poromaa, Department of Computing Science, */
00096 /*     Umea University, S-901 87 Umea, Sweden. */
00097 
00098 /*  ===================================================================== */
00099 
00100 /*     .. Parameters .. */
00101 /*     .. */
00102 /*     .. Local Scalars .. */
00103 /*     .. */
00104 /*     .. External Subroutines .. */
00105 /*     .. */
00106 /*     .. External Functions .. */
00107 /*     .. */
00108 /*     .. Intrinsic Functions .. */
00109 /*     .. */
00110 /*     .. Executable Statements .. */
00111 
00112 /*     Set constants to control overflow */
00113 
00114     /* Parameter adjustments */
00115     a_dim1 = *lda;
00116     a_offset = 1 + a_dim1;
00117     a -= a_offset;
00118     --ipiv;
00119     --jpiv;
00120 
00121     /* Function Body */
00122     *info = 0;
00123     eps = dlamch_("P");
00124     smlnum = dlamch_("S") / eps;
00125     bignum = 1. / smlnum;
00126     dlabad_(&smlnum, &bignum);
00127 
00128 /*     Factorize A using complete pivoting. */
00129 /*     Set pivots less than SMIN to SMIN. */
00130 
00131     i__1 = *n - 1;
00132     for (i__ = 1; i__ <= i__1; ++i__) {
00133 
00134 /*        Find max element in matrix A */
00135 
00136         xmax = 0.;
00137         i__2 = *n;
00138         for (ip = i__; ip <= i__2; ++ip) {
00139             i__3 = *n;
00140             for (jp = i__; jp <= i__3; ++jp) {
00141                 if ((d__1 = a[ip + jp * a_dim1], abs(d__1)) >= xmax) {
00142                     xmax = (d__1 = a[ip + jp * a_dim1], abs(d__1));
00143                     ipv = ip;
00144                     jpv = jp;
00145                 }
00146 /* L10: */
00147             }
00148 /* L20: */
00149         }
00150         if (i__ == 1) {
00151 /* Computing MAX */
00152             d__1 = eps * xmax;
00153             smin = max(d__1,smlnum);
00154         }
00155 
00156 /*        Swap rows */
00157 
00158         if (ipv != i__) {
00159             dswap_(n, &a[ipv + a_dim1], lda, &a[i__ + a_dim1], lda);
00160         }
00161         ipiv[i__] = ipv;
00162 
00163 /*        Swap columns */
00164 
00165         if (jpv != i__) {
00166             dswap_(n, &a[jpv * a_dim1 + 1], &c__1, &a[i__ * a_dim1 + 1], &
00167                     c__1);
00168         }
00169         jpiv[i__] = jpv;
00170 
00171 /*        Check for singularity */
00172 
00173         if ((d__1 = a[i__ + i__ * a_dim1], abs(d__1)) < smin) {
00174             *info = i__;
00175             a[i__ + i__ * a_dim1] = smin;
00176         }
00177         i__2 = *n;
00178         for (j = i__ + 1; j <= i__2; ++j) {
00179             a[j + i__ * a_dim1] /= a[i__ + i__ * a_dim1];
00180 /* L30: */
00181         }
00182         i__2 = *n - i__;
00183         i__3 = *n - i__;
00184         dger_(&i__2, &i__3, &c_b10, &a[i__ + 1 + i__ * a_dim1], &c__1, &a[i__ 
00185                 + (i__ + 1) * a_dim1], lda, &a[i__ + 1 + (i__ + 1) * a_dim1], 
00186                 lda);
00187 /* L40: */
00188     }
00189 
00190     if ((d__1 = a[*n + *n * a_dim1], abs(d__1)) < smin) {
00191         *info = *n;
00192         a[*n + *n * a_dim1] = smin;
00193     }
00194 
00195     return 0;
00196 
00197 /*     End of DGETC2 */
00198 
00199 } /* dgetc2_ */


swiftnav
Author(s):
autogenerated on Sat Jun 8 2019 18:55:45