annotate ext/clapack/src/dgetrf.c @ 209:ccd2019190bf msvc

Some MSVC fixes, including (temporarily, probably) renaming the FFT source file to avoid getting it mixed up with the Vamp SDK one in our object dir
author Chris Cannam
date Thu, 01 Feb 2018 16:34:08 +0000
parents 45330e0d2819
children
rev   line source
Chris@202 1 /* dgetrf.f -- translated by f2c (version 20061008).
Chris@202 2 You must link the resulting object file with libf2c:
Chris@202 3 on Microsoft Windows system, link with libf2c.lib;
Chris@202 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
Chris@202 5 or, if you install libf2c.a in a standard place, with -lf2c -lm
Chris@202 6 -- in that order, at the end of the command line, as in
Chris@202 7 cc *.o -lf2c -lm
Chris@202 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
Chris@202 9
Chris@202 10 http://www.netlib.org/f2c/libf2c.zip
Chris@202 11 */
Chris@202 12
Chris@202 13 #include "f2c.h"
Chris@202 14 #include "blaswrap.h"
Chris@202 15
Chris@202 16 /* Table of constant values */
Chris@202 17
Chris@202 18 static integer c__1 = 1;
Chris@202 19 static integer c_n1 = -1;
Chris@202 20 static doublereal c_b16 = 1.;
Chris@202 21 static doublereal c_b19 = -1.;
Chris@202 22
Chris@202 23 /* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer *
Chris@202 24 lda, integer *ipiv, integer *info)
Chris@202 25 {
Chris@202 26 /* System generated locals */
Chris@202 27 integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
Chris@202 28
Chris@202 29 /* Local variables */
Chris@202 30 integer i__, j, jb, nb;
Chris@202 31 extern /* Subroutine */ int dgemm_(char *, char *, integer *, integer *,
Chris@202 32 integer *, doublereal *, doublereal *, integer *, doublereal *,
Chris@202 33 integer *, doublereal *, doublereal *, integer *);
Chris@202 34 integer iinfo;
Chris@202 35 extern /* Subroutine */ int dtrsm_(char *, char *, char *, char *,
Chris@202 36 integer *, integer *, doublereal *, doublereal *, integer *,
Chris@202 37 doublereal *, integer *), dgetf2_(
Chris@202 38 integer *, integer *, doublereal *, integer *, integer *, integer
Chris@202 39 *), xerbla_(char *, integer *);
Chris@202 40 extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
Chris@202 41 integer *, integer *);
Chris@202 42 extern /* Subroutine */ int dlaswp_(integer *, doublereal *, integer *,
Chris@202 43 integer *, integer *, integer *, integer *);
Chris@202 44
Chris@202 45
Chris@202 46 /* -- LAPACK routine (version 3.2) -- */
Chris@202 47 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
Chris@202 48 /* November 2006 */
Chris@202 49
Chris@202 50 /* .. Scalar Arguments .. */
Chris@202 51 /* .. */
Chris@202 52 /* .. Array Arguments .. */
Chris@202 53 /* .. */
Chris@202 54
Chris@202 55 /* Purpose */
Chris@202 56 /* ======= */
Chris@202 57
Chris@202 58 /* DGETRF computes an LU factorization of a general M-by-N matrix A */
Chris@202 59 /* using partial pivoting with row interchanges. */
Chris@202 60
Chris@202 61 /* The factorization has the form */
Chris@202 62 /* A = P * L * U */
Chris@202 63 /* where P is a permutation matrix, L is lower triangular with unit */
Chris@202 64 /* diagonal elements (lower trapezoidal if m > n), and U is upper */
Chris@202 65 /* triangular (upper trapezoidal if m < n). */
Chris@202 66
Chris@202 67 /* This is the right-looking Level 3 BLAS version of the algorithm. */
Chris@202 68
Chris@202 69 /* Arguments */
Chris@202 70 /* ========= */
Chris@202 71
Chris@202 72 /* M (input) INTEGER */
Chris@202 73 /* The number of rows of the matrix A. M >= 0. */
Chris@202 74
Chris@202 75 /* N (input) INTEGER */
Chris@202 76 /* The number of columns of the matrix A. N >= 0. */
Chris@202 77
Chris@202 78 /* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
Chris@202 79 /* On entry, the M-by-N matrix to be factored. */
Chris@202 80 /* On exit, the factors L and U from the factorization */
Chris@202 81 /* A = P*L*U; the unit diagonal elements of L are not stored. */
Chris@202 82
Chris@202 83 /* LDA (input) INTEGER */
Chris@202 84 /* The leading dimension of the array A. LDA >= max(1,M). */
Chris@202 85
Chris@202 86 /* IPIV (output) INTEGER array, dimension (min(M,N)) */
Chris@202 87 /* The pivot indices; for 1 <= i <= min(M,N), row i of the */
Chris@202 88 /* matrix was interchanged with row IPIV(i). */
Chris@202 89
Chris@202 90 /* INFO (output) INTEGER */
Chris@202 91 /* = 0: successful exit */
Chris@202 92 /* < 0: if INFO = -i, the i-th argument had an illegal value */
Chris@202 93 /* > 0: if INFO = i, U(i,i) is exactly zero. The factorization */
Chris@202 94 /* has been completed, but the factor U is exactly */
Chris@202 95 /* singular, and division by zero will occur if it is used */
Chris@202 96 /* to solve a system of equations. */
Chris@202 97
Chris@202 98 /* ===================================================================== */
Chris@202 99
Chris@202 100 /* .. Parameters .. */
Chris@202 101 /* .. */
Chris@202 102 /* .. Local Scalars .. */
Chris@202 103 /* .. */
Chris@202 104 /* .. External Subroutines .. */
Chris@202 105 /* .. */
Chris@202 106 /* .. External Functions .. */
Chris@202 107 /* .. */
Chris@202 108 /* .. Intrinsic Functions .. */
Chris@202 109 /* .. */
Chris@202 110 /* .. Executable Statements .. */
Chris@202 111
Chris@202 112 /* Test the input parameters. */
Chris@202 113
Chris@202 114 /* Parameter adjustments */
Chris@202 115 a_dim1 = *lda;
Chris@202 116 a_offset = 1 + a_dim1;
Chris@202 117 a -= a_offset;
Chris@202 118 --ipiv;
Chris@202 119
Chris@202 120 /* Function Body */
Chris@202 121 *info = 0;
Chris@202 122 if (*m < 0) {
Chris@202 123 *info = -1;
Chris@202 124 } else if (*n < 0) {
Chris@202 125 *info = -2;
Chris@202 126 } else if (*lda < max(1,*m)) {
Chris@202 127 *info = -4;
Chris@202 128 }
Chris@202 129 if (*info != 0) {
Chris@202 130 i__1 = -(*info);
Chris@202 131 xerbla_("DGETRF", &i__1);
Chris@202 132 return 0;
Chris@202 133 }
Chris@202 134
Chris@202 135 /* Quick return if possible */
Chris@202 136
Chris@202 137 if (*m == 0 || *n == 0) {
Chris@202 138 return 0;
Chris@202 139 }
Chris@202 140
Chris@202 141 /* Determine the block size for this environment. */
Chris@202 142
Chris@202 143 nb = ilaenv_(&c__1, "DGETRF", " ", m, n, &c_n1, &c_n1);
Chris@202 144 if (nb <= 1 || nb >= min(*m,*n)) {
Chris@202 145
Chris@202 146 /* Use unblocked code. */
Chris@202 147
Chris@202 148 dgetf2_(m, n, &a[a_offset], lda, &ipiv[1], info);
Chris@202 149 } else {
Chris@202 150
Chris@202 151 /* Use blocked code. */
Chris@202 152
Chris@202 153 i__1 = min(*m,*n);
Chris@202 154 i__2 = nb;
Chris@202 155 for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
Chris@202 156 /* Computing MIN */
Chris@202 157 i__3 = min(*m,*n) - j + 1;
Chris@202 158 jb = min(i__3,nb);
Chris@202 159
Chris@202 160 /* Factor diagonal and subdiagonal blocks and test for exact */
Chris@202 161 /* singularity. */
Chris@202 162
Chris@202 163 i__3 = *m - j + 1;
Chris@202 164 dgetf2_(&i__3, &jb, &a[j + j * a_dim1], lda, &ipiv[j], &iinfo);
Chris@202 165
Chris@202 166 /* Adjust INFO and the pivot indices. */
Chris@202 167
Chris@202 168 if (*info == 0 && iinfo > 0) {
Chris@202 169 *info = iinfo + j - 1;
Chris@202 170 }
Chris@202 171 /* Computing MIN */
Chris@202 172 i__4 = *m, i__5 = j + jb - 1;
Chris@202 173 i__3 = min(i__4,i__5);
Chris@202 174 for (i__ = j; i__ <= i__3; ++i__) {
Chris@202 175 ipiv[i__] = j - 1 + ipiv[i__];
Chris@202 176 /* L10: */
Chris@202 177 }
Chris@202 178
Chris@202 179 /* Apply interchanges to columns 1:J-1. */
Chris@202 180
Chris@202 181 i__3 = j - 1;
Chris@202 182 i__4 = j + jb - 1;
Chris@202 183 dlaswp_(&i__3, &a[a_offset], lda, &j, &i__4, &ipiv[1], &c__1);
Chris@202 184
Chris@202 185 if (j + jb <= *n) {
Chris@202 186
Chris@202 187 /* Apply interchanges to columns J+JB:N. */
Chris@202 188
Chris@202 189 i__3 = *n - j - jb + 1;
Chris@202 190 i__4 = j + jb - 1;
Chris@202 191 dlaswp_(&i__3, &a[(j + jb) * a_dim1 + 1], lda, &j, &i__4, &
Chris@202 192 ipiv[1], &c__1);
Chris@202 193
Chris@202 194 /* Compute block row of U. */
Chris@202 195
Chris@202 196 i__3 = *n - j - jb + 1;
Chris@202 197 dtrsm_("Left", "Lower", "No transpose", "Unit", &jb, &i__3, &
Chris@202 198 c_b16, &a[j + j * a_dim1], lda, &a[j + (j + jb) *
Chris@202 199 a_dim1], lda);
Chris@202 200 if (j + jb <= *m) {
Chris@202 201
Chris@202 202 /* Update trailing submatrix. */
Chris@202 203
Chris@202 204 i__3 = *m - j - jb + 1;
Chris@202 205 i__4 = *n - j - jb + 1;
Chris@202 206 dgemm_("No transpose", "No transpose", &i__3, &i__4, &jb,
Chris@202 207 &c_b19, &a[j + jb + j * a_dim1], lda, &a[j + (j +
Chris@202 208 jb) * a_dim1], lda, &c_b16, &a[j + jb + (j + jb) *
Chris@202 209 a_dim1], lda);
Chris@202 210 }
Chris@202 211 }
Chris@202 212 /* L20: */
Chris@202 213 }
Chris@202 214 }
Chris@202 215 return 0;
Chris@202 216
Chris@202 217 /* End of DGETRF */
Chris@202 218
Chris@202 219 } /* dgetrf_ */