git.maemo.org Git - opencv/blob - 3rdparty/lapack/sgetrf.c

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6 static integer c_n1 = -1;
   7 static real c_b16 = 1.f;
   8 static real c_b19 = -1.f;
   9
  10 /* Subroutine */ int sgetrf_(integer *m, integer *n, real *a, integer *lda,
  11         integer *ipiv, integer *info)
  12 {
  13     /* System generated locals */
  14     integer a_dim1, a_offset, i__1, i__2, i__3, i__4, i__5;
  15
  16     /* Local variables */
  17     integer i__, j, jb, nb, iinfo;
  18     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
  19             integer *, real *, real *, integer *, real *, integer *, real *,
  20             real *, integer *), strsm_(char *, char *, char *,
  21              char *, integer *, integer *, real *, real *, integer *, real *,
  22             integer *), sgetf2_(integer *,
  23             integer *, real *, integer *, integer *, integer *), xerbla_(char
  24             *, integer *);
  25     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  26             integer *, integer *);
  27     extern /* Subroutine */ int slaswp_(integer *, real *, integer *, integer
  28             *, integer *, integer *, integer *);
  29
  30
  31 /*  -- LAPACK routine (version 3.1) -- */
  32 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  33 /*     November 2006 */
  34
  35 /*     .. Scalar Arguments .. */
  36 /*     .. */
  37 /*     .. Array Arguments .. */
  38 /*     .. */
  39
  40 /*  Purpose */
  41 /*  ======= */
  42
  43 /*  SGETRF computes an LU factorization of a general M-by-N matrix A */
  44 /*  using partial pivoting with row interchanges. */
  45
  46 /*  The factorization has the form */
  47 /*     A = P * L * U */
  48 /*  where P is a permutation matrix, L is lower triangular with unit */
  49 /*  diagonal elements (lower trapezoidal if m > n), and U is upper */
  50 /*  triangular (upper trapezoidal if m < n). */
  51
  52 /*  This is the right-looking Level 3 BLAS version of the algorithm. */
  53
  54 /*  Arguments */
  55 /*  ========= */
  56
  57 /*  M       (input) INTEGER */
  58 /*          The number of rows of the matrix A.  M >= 0. */
  59
  60 /*  N       (input) INTEGER */
  61 /*          The number of columns of the matrix A.  N >= 0. */
  62
  63 /*  A       (input/output) REAL array, dimension (LDA,N) */
  64 /*          On entry, the M-by-N matrix to be factored. */
  65 /*          On exit, the factors L and U from the factorization */
  66 /*          A = P*L*U; the unit diagonal elements of L are not stored. */
  67
  68 /*  LDA     (input) INTEGER */
  69 /*          The leading dimension of the array A.  LDA >= max(1,M). */
  70
  71 /*  IPIV    (output) INTEGER array, dimension (min(M,N)) */
  72 /*          The pivot indices; for 1 <= i <= min(M,N), row i of the */
  73 /*          matrix was interchanged with row IPIV(i). */
  74
  75 /*  INFO    (output) INTEGER */
  76 /*          = 0:  successful exit */
  77 /*          < 0:  if INFO = -i, the i-th argument had an illegal value */
  78 /*          > 0:  if INFO = i, U(i,i) is exactly zero. The factorization */
  79 /*                has been completed, but the factor U is exactly */
  80 /*                singular, and division by zero will occur if it is used */
  81 /*                to solve a system of equations. */
  82
  83 /*  ===================================================================== */
  84
  85 /*     .. Parameters .. */
  86 /*     .. */
  87 /*     .. Local Scalars .. */
  88 /*     .. */
  89 /*     .. External Subroutines .. */
  90 /*     .. */
  91 /*     .. External Functions .. */
  92 /*     .. */
  93 /*     .. Intrinsic Functions .. */
  94 /*     .. */
  95 /*     .. Executable Statements .. */
  96
  97 /*     Test the input parameters. */
  98
  99     /* Parameter adjustments */
 100     a_dim1 = *lda;
 101     a_offset = 1 + a_dim1;
 102     a -= a_offset;
 103     --ipiv;
 104
 105     /* Function Body */
 106     *info = 0;
 107     if (*m < 0) {
 108         *info = -1;
 109     } else if (*n < 0) {
 110         *info = -2;
 111     } else if (*lda < max(1,*m)) {
 112         *info = -4;
 113     }
 114     if (*info != 0) {
 115         i__1 = -(*info);
 116         xerbla_("SGETRF", &i__1);
 117         return 0;
 118     }
 119
 120 /*     Quick return if possible */
 121
 122     if (*m == 0 || *n == 0) {
 123         return 0;
 124     }
 125
 126 /*     Determine the block size for this environment. */
 127
 128     nb = ilaenv_(&c__1, "SGETRF", " ", m, n, &c_n1, &c_n1);
 129     if (nb <= 1 || nb >= min(*m,*n)) {
 130
 131 /*        Use unblocked code. */
 132
 133         sgetf2_(m, n, &a[a_offset], lda, &ipiv[1], info);
 134     } else {
 135
 136 /*        Use blocked code. */
 137
 138         i__1 = min(*m,*n);
 139         i__2 = nb;
 140         for (j = 1; i__2 < 0 ? j >= i__1 : j <= i__1; j += i__2) {
 141 /* Computing MIN */
 142             i__3 = min(*m,*n) - j + 1;
 143             jb = min(i__3,nb);
 144
 145 /*           Factor diagonal and subdiagonal blocks and test for exact */
 146 /*           singularity. */
 147
 148             i__3 = *m - j + 1;
 149             sgetf2_(&i__3, &jb, &a[j + j * a_dim1], lda, &ipiv[j], &iinfo);
 150
 151 /*           Adjust INFO and the pivot indices. */
 152
 153             if (*info == 0 && iinfo > 0) {
 154                 *info = iinfo + j - 1;
 155             }
 156 /* Computing MIN */
 157             i__4 = *m, i__5 = j + jb - 1;
 158             i__3 = min(i__4,i__5);
 159             for (i__ = j; i__ <= i__3; ++i__) {
 160                 ipiv[i__] = j - 1 + ipiv[i__];
 161 /* L10: */
 162             }
 163
 164 /*           Apply interchanges to columns 1:J-1. */
 165
 166             i__3 = j - 1;
 167             i__4 = j + jb - 1;
 168             slaswp_(&i__3, &a[a_offset], lda, &j, &i__4, &ipiv[1], &c__1);
 169
 170             if (j + jb <= *n) {
 171
 172 /*              Apply interchanges to columns J+JB:N. */
 173
 174                 i__3 = *n - j - jb + 1;
 175                 i__4 = j + jb - 1;
 176                 slaswp_(&i__3, &a[(j + jb) * a_dim1 + 1], lda, &j, &i__4, &
 177                         ipiv[1], &c__1);
 178
 179 /*              Compute block row of U. */
 180
 181                 i__3 = *n - j - jb + 1;
 182                 strsm_("Left", "Lower", "No transpose", "Unit", &jb, &i__3, &
 183                         c_b16, &a[j + j * a_dim1], lda, &a[j + (j + jb) *
 184                         a_dim1], lda);
 185                 if (j + jb <= *m) {
 186
 187 /*                 Update trailing submatrix. */
 188
 189                     i__3 = *m - j - jb + 1;
 190                     i__4 = *n - j - jb + 1;
 191                     sgemm_("No transpose", "No transpose", &i__3, &i__4, &jb,
 192                             &c_b19, &a[j + jb + j * a_dim1], lda, &a[j + (j +
 193                             jb) * a_dim1], lda, &c_b16, &a[j + jb + (j + jb) *
 194                              a_dim1], lda);
 195                 }
 196             }
 197 /* L20: */
 198         }
 199     }
 200     return 0;
 201
 202 /*     End of SGETRF */
 203
 204 } /* sgetrf_ */