git.maemo.org Git - opencv/blob - 3rdparty/lapack/slauum.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_b15 = 1.f;
   8
   9 /* Subroutine */ int slauum_(char *uplo, integer *n, real *a, integer *lda,
  10         integer *info)
  11 {
  12     /* System generated locals */
  13     integer a_dim1, a_offset, i__1, i__2, i__3, i__4;
  14
  15     /* Local variables */
  16     integer i__, ib, nb;
  17     extern logical lsame_(char *, char *);
  18     extern /* Subroutine */ int sgemm_(char *, char *, integer *, integer *,
  19             integer *, real *, real *, integer *, real *, integer *, real *,
  20             real *, integer *);
  21     logical upper;
  22     extern /* Subroutine */ int strmm_(char *, char *, char *, char *,
  23             integer *, integer *, real *, real *, integer *, real *, integer *
  24 ), ssyrk_(char *, char *, integer
  25             *, integer *, real *, real *, integer *, real *, real *, integer *
  26 ), slauu2_(char *, integer *, real *, integer *,
  27             integer *), xerbla_(char *, integer *);
  28     extern integer ilaenv_(integer *, char *, char *, integer *, integer *,
  29             integer *, integer *);
  30
  31
  32 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  33 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  34 /*     November 2006 */
  35
  36 /*     .. Scalar Arguments .. */
  37 /*     .. */
  38 /*     .. Array Arguments .. */
  39 /*     .. */
  40
  41 /*  Purpose */
  42 /*  ======= */
  43
  44 /*  SLAUUM computes the product U * U' or L' * L, where the triangular */
  45 /*  factor U or L is stored in the upper or lower triangular part of */
  46 /*  the array A. */
  47
  48 /*  If UPLO = 'U' or 'u' then the upper triangle of the result is stored, */
  49 /*  overwriting the factor U in A. */
  50 /*  If UPLO = 'L' or 'l' then the lower triangle of the result is stored, */
  51 /*  overwriting the factor L in A. */
  52
  53 /*  This is the blocked form of the algorithm, calling Level 3 BLAS. */
  54
  55 /*  Arguments */
  56 /*  ========= */
  57
  58 /*  UPLO    (input) CHARACTER*1 */
  59 /*          Specifies whether the triangular factor stored in the array A */
  60 /*          is upper or lower triangular: */
  61 /*          = 'U':  Upper triangular */
  62 /*          = 'L':  Lower triangular */
  63
  64 /*  N       (input) INTEGER */
  65 /*          The order of the triangular factor U or L.  N >= 0. */
  66
  67 /*  A       (input/output) REAL array, dimension (LDA,N) */
  68 /*          On entry, the triangular factor U or L. */
  69 /*          On exit, if UPLO = 'U', the upper triangle of A is */
  70 /*          overwritten with the upper triangle of the product U * U'; */
  71 /*          if UPLO = 'L', the lower triangle of A is overwritten with */
  72 /*          the lower triangle of the product L' * L. */
  73
  74 /*  LDA     (input) INTEGER */
  75 /*          The leading dimension of the array A.  LDA >= max(1,N). */
  76
  77 /*  INFO    (output) INTEGER */
  78 /*          = 0: successful exit */
  79 /*          < 0: if INFO = -k, the k-th argument had an illegal value */
  80
  81 /*  ===================================================================== */
  82
  83 /*     .. Parameters .. */
  84 /*     .. */
  85 /*     .. Local Scalars .. */
  86 /*     .. */
  87 /*     .. External Functions .. */
  88 /*     .. */
  89 /*     .. External Subroutines .. */
  90 /*     .. */
  91 /*     .. Intrinsic Functions .. */
  92 /*     .. */
  93 /*     .. Executable Statements .. */
  94
  95 /*     Test the input parameters. */
  96
  97     /* Parameter adjustments */
  98     a_dim1 = *lda;
  99     a_offset = 1 + a_dim1;
 100     a -= a_offset;
 101
 102     /* Function Body */
 103     *info = 0;
 104     upper = lsame_(uplo, "U");
 105     if (! upper && ! lsame_(uplo, "L")) {
 106         *info = -1;
 107     } else if (*n < 0) {
 108         *info = -2;
 109     } else if (*lda < max(1,*n)) {
 110         *info = -4;
 111     }
 112     if (*info != 0) {
 113         i__1 = -(*info);
 114         xerbla_("SLAUUM", &i__1);
 115         return 0;
 116     }
 117
 118 /*     Quick return if possible */
 119
 120     if (*n == 0) {
 121         return 0;
 122     }
 123
 124 /*     Determine the block size for this environment. */
 125
 126     nb = ilaenv_(&c__1, "SLAUUM", uplo, n, &c_n1, &c_n1, &c_n1);
 127
 128     if (nb <= 1 || nb >= *n) {
 129
 130 /*        Use unblocked code */
 131
 132         slauu2_(uplo, n, &a[a_offset], lda, info);
 133     } else {
 134
 135 /*        Use blocked code */
 136
 137         if (upper) {
 138
 139 /*           Compute the product U * U'. */
 140
 141             i__1 = *n;
 142             i__2 = nb;
 143             for (i__ = 1; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) {
 144 /* Computing MIN */
 145                 i__3 = nb, i__4 = *n - i__ + 1;
 146                 ib = min(i__3,i__4);
 147                 i__3 = i__ - 1;
 148                 strmm_("Right", "Upper", "Transpose", "Non-unit", &i__3, &ib,
 149                         &c_b15, &a[i__ + i__ * a_dim1], lda, &a[i__ * a_dim1
 150                         + 1], lda)
 151                         ;
 152                 slauu2_("Upper", &ib, &a[i__ + i__ * a_dim1], lda, info);
 153                 if (i__ + ib <= *n) {
 154                     i__3 = i__ - 1;
 155                     i__4 = *n - i__ - ib + 1;
 156                     sgemm_("No transpose", "Transpose", &i__3, &ib, &i__4, &
 157                             c_b15, &a[(i__ + ib) * a_dim1 + 1], lda, &a[i__ +
 158                             (i__ + ib) * a_dim1], lda, &c_b15, &a[i__ *
 159                             a_dim1 + 1], lda);
 160                     i__3 = *n - i__ - ib + 1;
 161                     ssyrk_("Upper", "No transpose", &ib, &i__3, &c_b15, &a[
 162                             i__ + (i__ + ib) * a_dim1], lda, &c_b15, &a[i__ +
 163                             i__ * a_dim1], lda);
 164                 }
 165 /* L10: */
 166             }
 167         } else {
 168
 169 /*           Compute the product L' * L. */
 170
 171             i__2 = *n;
 172             i__1 = nb;
 173             for (i__ = 1; i__1 < 0 ? i__ >= i__2 : i__ <= i__2; i__ += i__1) {
 174 /* Computing MIN */
 175                 i__3 = nb, i__4 = *n - i__ + 1;
 176                 ib = min(i__3,i__4);
 177                 i__3 = i__ - 1;
 178                 strmm_("Left", "Lower", "Transpose", "Non-unit", &ib, &i__3, &
 179                         c_b15, &a[i__ + i__ * a_dim1], lda, &a[i__ + a_dim1],
 180                         lda);
 181                 slauu2_("Lower", &ib, &a[i__ + i__ * a_dim1], lda, info);
 182                 if (i__ + ib <= *n) {
 183                     i__3 = i__ - 1;
 184                     i__4 = *n - i__ - ib + 1;
 185                     sgemm_("Transpose", "No transpose", &ib, &i__3, &i__4, &
 186                             c_b15, &a[i__ + ib + i__ * a_dim1], lda, &a[i__ +
 187                             ib + a_dim1], lda, &c_b15, &a[i__ + a_dim1], lda);
 188                     i__3 = *n - i__ - ib + 1;
 189                     ssyrk_("Lower", "Transpose", &ib, &i__3, &c_b15, &a[i__ +
 190                             ib + i__ * a_dim1], lda, &c_b15, &a[i__ + i__ *
 191                             a_dim1], lda);
 192                 }
 193 /* L20: */
 194             }
 195         }
 196     }
 197
 198     return 0;
 199
 200 /*     End of SLAUUM */
 201
 202 } /* slauum_ */