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

   1 #include "clapack.h"
   2
   3 /* Table of constant values */
   4
   5 static integer c__1 = 1;
   6
   7 doublereal slansy_(char *norm, char *uplo, integer *n, real *a, integer *lda,
   8         real *work)
   9 {
  10     /* System generated locals */
  11     integer a_dim1, a_offset, i__1, i__2;
  12     real ret_val, r__1, r__2, r__3;
  13
  14     /* Builtin functions */
  15     double sqrt(doublereal);
  16
  17     /* Local variables */
  18     integer i__, j;
  19     real sum, absa, scale;
  20     extern logical lsame_(char *, char *);
  21     real value;
  22     extern /* Subroutine */ int slassq_(integer *, real *, integer *, real *,
  23             real *);
  24
  25
  26 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  27 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  28 /*     November 2006 */
  29
  30 /*     .. Scalar Arguments .. */
  31 /*     .. */
  32 /*     .. Array Arguments .. */
  33 /*     .. */
  34
  35 /*  Purpose */
  36 /*  ======= */
  37
  38 /*  SLANSY  returns the value of the one norm,  or the Frobenius norm, or */
  39 /*  the  infinity norm,  or the  element of  largest absolute value  of a */
  40 /*  real symmetric matrix A. */
  41
  42 /*  Description */
  43 /*  =========== */
  44
  45 /*  SLANSY returns the value */
  46
  47 /*     SLANSY = ( max(abs(A(i,j))), NORM = 'M' or 'm' */
  48 /*              ( */
  49 /*              ( norm1(A),         NORM = '1', 'O' or 'o' */
  50 /*              ( */
  51 /*              ( normI(A),         NORM = 'I' or 'i' */
  52 /*              ( */
  53 /*              ( normF(A),         NORM = 'F', 'f', 'E' or 'e' */
  54
  55 /*  where  norm1  denotes the  one norm of a matrix (maximum column sum), */
  56 /*  normI  denotes the  infinity norm  of a matrix  (maximum row sum) and */
  57 /*  normF  denotes the  Frobenius norm of a matrix (square root of sum of */
  58 /*  squares).  Note that  max(abs(A(i,j)))  is not a consistent matrix norm. */
  59
  60 /*  Arguments */
  61 /*  ========= */
  62
  63 /*  NORM    (input) CHARACTER*1 */
  64 /*          Specifies the value to be returned in SLANSY as described */
  65 /*          above. */
  66
  67 /*  UPLO    (input) CHARACTER*1 */
  68 /*          Specifies whether the upper or lower triangular part of the */
  69 /*          symmetric matrix A is to be referenced. */
  70 /*          = 'U':  Upper triangular part of A is referenced */
  71 /*          = 'L':  Lower triangular part of A is referenced */
  72
  73 /*  N       (input) INTEGER */
  74 /*          The order of the matrix A.  N >= 0.  When N = 0, SLANSY is */
  75 /*          set to zero. */
  76
  77 /*  A       (input) REAL array, dimension (LDA,N) */
  78 /*          The symmetric matrix A.  If UPLO = 'U', the leading n by n */
  79 /*          upper triangular part of A contains the upper triangular part */
  80 /*          of the matrix A, and the strictly lower triangular part of A */
  81 /*          is not referenced.  If UPLO = 'L', the leading n by n lower */
  82 /*          triangular part of A contains the lower triangular part of */
  83 /*          the matrix A, and the strictly upper triangular part of A is */
  84 /*          not referenced. */
  85
  86 /*  LDA     (input) INTEGER */
  87 /*          The leading dimension of the array A.  LDA >= max(N,1). */
  88
  89 /*  WORK    (workspace) REAL array, dimension (MAX(1,LWORK)), */
  90 /*          where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise, */
  91 /*          WORK is not referenced. */
  92
  93 /* ===================================================================== */
  94
  95 /*     .. Parameters .. */
  96 /*     .. */
  97 /*     .. Local Scalars .. */
  98 /*     .. */
  99 /*     .. External Subroutines .. */
 100 /*     .. */
 101 /*     .. External Functions .. */
 102 /*     .. */
 103 /*     .. Intrinsic Functions .. */
 104 /*     .. */
 105 /*     .. Executable Statements .. */
 106
 107     /* Parameter adjustments */
 108     a_dim1 = *lda;
 109     a_offset = 1 + a_dim1;
 110     a -= a_offset;
 111     --work;
 112
 113     /* Function Body */
 114     if (*n == 0) {
 115         value = 0.f;
 116     } else if (lsame_(norm, "M")) {
 117
 118 /*        Find max(abs(A(i,j))). */
 119
 120         value = 0.f;
 121         if (lsame_(uplo, "U")) {
 122             i__1 = *n;
 123             for (j = 1; j <= i__1; ++j) {
 124                 i__2 = j;
 125                 for (i__ = 1; i__ <= i__2; ++i__) {
 126 /* Computing MAX */
 127                     r__2 = value, r__3 = (r__1 = a[i__ + j * a_dim1], dabs(
 128                             r__1));
 129                     value = dmax(r__2,r__3);
 130 /* L10: */
 131                 }
 132 /* L20: */
 133             }
 134         } else {
 135             i__1 = *n;
 136             for (j = 1; j <= i__1; ++j) {
 137                 i__2 = *n;
 138                 for (i__ = j; i__ <= i__2; ++i__) {
 139 /* Computing MAX */
 140                     r__2 = value, r__3 = (r__1 = a[i__ + j * a_dim1], dabs(
 141                             r__1));
 142                     value = dmax(r__2,r__3);
 143 /* L30: */
 144                 }
 145 /* L40: */
 146             }
 147         }
 148     } else if (lsame_(norm, "I") || lsame_(norm, "O") || *(unsigned char *)norm == '1') {
 149
 150 /*        Find normI(A) ( = norm1(A), since A is symmetric). */
 151
 152         value = 0.f;
 153         if (lsame_(uplo, "U")) {
 154             i__1 = *n;
 155             for (j = 1; j <= i__1; ++j) {
 156                 sum = 0.f;
 157                 i__2 = j - 1;
 158                 for (i__ = 1; i__ <= i__2; ++i__) {
 159                     absa = (r__1 = a[i__ + j * a_dim1], dabs(r__1));
 160                     sum += absa;
 161                     work[i__] += absa;
 162 /* L50: */
 163                 }
 164                 work[j] = sum + (r__1 = a[j + j * a_dim1], dabs(r__1));
 165 /* L60: */
 166             }
 167             i__1 = *n;
 168             for (i__ = 1; i__ <= i__1; ++i__) {
 169 /* Computing MAX */
 170                 r__1 = value, r__2 = work[i__];
 171                 value = dmax(r__1,r__2);
 172 /* L70: */
 173             }
 174         } else {
 175             i__1 = *n;
 176             for (i__ = 1; i__ <= i__1; ++i__) {
 177                 work[i__] = 0.f;
 178 /* L80: */
 179             }
 180             i__1 = *n;
 181             for (j = 1; j <= i__1; ++j) {
 182                 sum = work[j] + (r__1 = a[j + j * a_dim1], dabs(r__1));
 183                 i__2 = *n;
 184                 for (i__ = j + 1; i__ <= i__2; ++i__) {
 185                     absa = (r__1 = a[i__ + j * a_dim1], dabs(r__1));
 186                     sum += absa;
 187                     work[i__] += absa;
 188 /* L90: */
 189                 }
 190                 value = dmax(value,sum);
 191 /* L100: */
 192             }
 193         }
 194     } else if (lsame_(norm, "F") || lsame_(norm, "E")) {
 195
 196 /*        Find normF(A). */
 197
 198         scale = 0.f;
 199         sum = 1.f;
 200         if (lsame_(uplo, "U")) {
 201             i__1 = *n;
 202             for (j = 2; j <= i__1; ++j) {
 203                 i__2 = j - 1;
 204                 slassq_(&i__2, &a[j * a_dim1 + 1], &c__1, &scale, &sum);
 205 /* L110: */
 206             }
 207         } else {
 208             i__1 = *n - 1;
 209             for (j = 1; j <= i__1; ++j) {
 210                 i__2 = *n - j;
 211                 slassq_(&i__2, &a[j + 1 + j * a_dim1], &c__1, &scale, &sum);
 212 /* L120: */
 213             }
 214         }
 215         sum *= 2;
 216         i__1 = *lda + 1;
 217         slassq_(n, &a[a_offset], &i__1, &scale, &sum);
 218         value = scale * sqrt(sum);
 219     }
 220
 221     ret_val = value;
 222     return ret_val;
 223
 224 /*     End of SLANSY */
 225
 226 } /* slansy_ */